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occur with the application of the nudging method to incorporate sea surface height 
observations (with associated vertical corrections to temperature and salinity) into a global 
eddy resolving ocean model. The first experiment used a previous model run as the 
observational field to determine if the assimilation technique, nudging, produced significant 
changes in the simulated fields to adjust the model to the observed fields when starting at a 
statistically different initial condition. The twin experiment has shown that the model does 
respond to the inclusion of the observed fields. Both the surface fields and subsurface 
fields have been adjusted towards these synthetic observations. The second experiment 
involved the use of a combined altimetric sea surface height anomaly field from the ERS-1 
and the T/P satellites. The surface height fields are extended vertically by using the Levitus 
94 monthly climatological fields. This dissertation has shown that assimilation of surface 
height data and an associated vertical adjustment to temperature and salinity, modifies both 
the surface and subsurface fields. Changes can be seen in both prognostic and diagnostic 
quantities (such as heat content and meridional overturning) while remaining dynamically 
consistent with the numerics of the model itself. Comparison of the simulated fields with 
in situ observations of temperature and salinity show that the model has adjusted towards 
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I. INTRODUCTION 



A. PURPOSE 

High resolution global ocean models now produce very realistic simulations in 
which currents are adequately represented in strength and width, eddies are generated in the 
correct ocean regions where they are seen in hydrographic data, and seasonal to interannual 
changes of the near surface flow compare favorably with altimetric measurements of sea 
surface height fields. These models simulate the narrow widths and strength of the global 
currents which transport the majority of heat and salt within and between basins. It is 
important that the currents be accurately represented if any reasonable estimate of the global 
climate changes are to be made. However, these models are deficient in several ways: 
timing of turbulent (eddy) activity does not always coincide with reality, the surface heat 
flux does not get transferred to subsurface layers in a realistic manner, and initial 
conditions determined from climatological data sets are not fully correct and most likely, 
need to be adjusted. Without resorting to modifying the numerics of the basic model, some 
of these deficiencies can be identified and overcome by assimilating surface and subsurface 
observational data into the model. Excellent global altimetric data sets from both ERS-1/2 
and TOPEX/POSEIDON (T/P) provide a necessary global observational data to complete 
the task. 

This thesis describes such an assimilation process and the physical and numerical 
results. For the first time, the assimilation process has been applied to a global, high 
resolution primitive equation model ocean model with a free surface. In addition, the 
assimilation technique of nudging has been extended to consider the variance of the 
observational data, along with adjusting both the temperature and salinity in the subsurface 
layers of the model. The analysis of the simulations discusses how assimilation of 
observational data can be used to identify weaknesses in the numerical model formulation 
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and suggests how the models can be improved. The analysis is both qualitatitve and 
quantitative, evaluating both point measurements and integrated quantities. 

B. OCEAN GENERAL CIRCULATION MODEL 

The model used in this research is the Semmer/Chervin global ocean model, 
referred to also as the Parallel Ocean Climate Model (POCM). It is a primitive equation 
(PE) model based on the early formulation of Bryan (1969), updated by Semtner (1974) 
and modified by Semtner and Chervin (1988, 1992) to run on a parallel vector computer 
with multiple processors. The model has 20 layers, uses a B-grid (Mesinger and Arakawa, 
1976), and a Mercator grid of a 1/4° average cell size. The model is forced by 3-day 
average wind stresses, computed and interpolated from twice daily European Centre for 
Medium Weather Forecasting (ECMWF) 10-meter winds. Climatological monthly heat 
fluxes (Barnier et al. 1994), interpolated to 3-day fields, are used and the salt fluxes are 
handled by relaxing to the Levitus 1994 surface salinity values {Levitus and Boyer, 1994). 
Subsurface restoring to observations occurs only in the vicinity of Gibraltar and in the 
high-latitude belts 58-65°N and 68-75°S. The model has a prognostic free surface 
(Killworth et al. 1991) to facilitate the nudging of the sea surface height (SSH) field 
towards the altimeter data and allow for realistic coastlines and bottom topography. Initial 
conditions are based on a 32-year 1/2° spinup from gridded data, equilibrated at 1/4° using 
1985-89 monthly 1000-mb winds. Detailed analysis of the accuracy of the large scale 
circulation simulated by POCM_4B (the name comes from it being the second run of the 
Semmer/Chervin POCM model at 1/4° average resolution) can be found in Stammer et al. 
(1996) and of the mesoscale in McClean etal. (1997). From the analysis of the nine- year 
run (1987-1995), it was found that the model is producing realistic circulation patterns on 
the large scale, but is somewhat deficient in eddy kinetic energy and the uptake of heat in 
the top layers. Tokmakian (1996) conducted a study to examine how well POCM 
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represented the time series of SSH rather than the statistical variability of SSH. Correlation 
values computed between the models and tide gauge data for the period of 1988 through 
1994 using the POCM have an average correlation over all 94 stations of 0.46 with a 
median value 0.50. Figure 1 shows several time series of model SSH and measurements 
made by tide gauges. In all three cases the model is reproducing the observed signal fairly 
well. The mismatch between model SSH and in situ data is partly the result of not having 
all the correct forcing included. In other words, if SSH includes the incorporation of an 
event forced at a time outside of the model run, then these events are not reflected in the 
model's estimate of SSH. A frequency analysis shows that the model is clearly 
representing the seasonal cycle in phase and amplitude accurately. The models' SSH 
anomalies represent, fairly well, the ocean's response to both local and remote wind driven 
events in a large part of the model. The model has active thermodynamics, allowing for 
temporal changes in surface heat and freshwater fluxes. 

The advantages of using a free-surface primitive equation model over one which 
has a rigid lid formulation can be seen by examining the basic formulation of the model and 
the corresponding prognostic variables of the model. In the formulation of a free surface 
model, the free surface, ii, is a prognostic variable. To determine the influence this 
prognostic variable has over the related diagnostic stream function of a rigid lid model in 
the assimilation of data, it needs to be determined where and to what extent the 
ageostrophic portion of the sea surface height contributes to the flow, both in the 
barotropic component of velocity and in the vertical velocity. 

The total u component of velocity can be written as u = Ug+ ua + Ue, where g is the 

geostrophic component, a is the ageostrophic component and e, the Ekman, or wind- 

driven, component. The free surface can be used to determine the geostrophic velocity, Ug 
g3ri 

= . This diagnostic variable is used to assimilate data into in a rigid lid model. With a 

free surface model, the field that is adjusted is the true SSH (t|), which is then coupled to 
the layers below by three equations: the continuity equation: 
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cm cm cm 




Figure 1 . Tide gauge measurement of sea level as compared to the P0CM-4B simulated SSH 
at a) Naos Island in the Gulf of Panama, b) at Townsville on the northeast coast of Australia, 
and c) at Christmas Island in the tropical Pacific. An annual harmonic has been removed 
from both signals. Noted on the plots axe the El Nino and Rossby wave events. The red line 
are the model estimates and the black line are the tide gauge measurements. 
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( 1 . 1 ) 



3u 8v 8w _ - 
3x 3y 3z ~ 

and the two momentum equations: 



and 



3u 3u 3u ^ 9 t| 



3v 3v 3v . 3 t| 



( 1 . 2 ) 



(1.3) 



By using a rigid lid formulation in a steady state, the first three terms of equations 

1 .2 and 1 .3 are assumed to be zero and the geostrophic definition of velocity results. If a 

steady state is assumed, then the relative importance of the 2nd and 3rd terms with respect 
to the fourth term needs to be evaluated using the model fields. These terms, u^ + v^ 
3v dv 

and u^ + v^ , can be computed using the model's top layer velocity field, with the 

Ekman component^ removed and then compared to the right hand side of the equation. 
The contribution of each of these terms to the total velocity can be estimated by looking at 
the ratio of the terms to the geostrophic term, fv or fu, respectively. In the tropics (± 10°), 
the ratio is -2.68 -646.51 for the u equation and 0.1 1120.45 for the v equation, suggesting 
that in both the u and v equations the non-linear terms dominate the geostrophic term. 
Outside of tropics, the terms have means of about zero, but with standard deviations of 
0.15 and 0.20, respectively, for the u and v equations. While over the whole area, the 
non-linear effects do not dominate the flow, there are some areas where they do affect the 
dynamics up to 20% outside the tropics. The global ratio for the u component is shown in 
Figure 2a, for values greater than 0.1 (or 10%) and for the v component in Figure 2b. 

In areas outside the tropics, such as in the region of the Gulf Stream and the 
Agulhas Current, the non-linear terms contribute to the total flow, especially at small 



^The non-Ekman flow in the top layer of the model can be determined by subtracting the Ekman flow, 



calculated from the wind stress: 
in the top layer. 



1 dXy 



(Xy is the y component of wind stress) from the total velocity 
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scales. This calls into question the general notion {Oschlies and Willebrand, 1996) that 
there is not much difference between assimilating SSH into a free surface model and into a 
rigid lid model. Thus it is physically better to use a free surface model to assimilate 
changes in SSH (t|) , rather than using a rigid lid model, so that the non-linear adjustments 
to T| are not ignored. 

C. ASSIMILATION METHODS 

This thesis concerns the assimilation of data into the above global model; therefore, 
some background to ocean assimilation is given. Large amounts of time and effort have 
been expended to investigate techniques for inserting altimeter data into ocean models. 
Ghil and Malanotte-Rizzoli (1991) and Malanotte-Rizzoli (1996) provide a general review 
of the current state of the assimilation using ocean models. Assimilation of data into 
models has been performed on a variety of models, quasi-geostrophic (QG), PE (layer and 
with active thermodynamics) and isopycnal models. Basically, there are two types of 
assimilation methods, sequential and variational. Sequential methods range from the very 
simple method of nudging to the much more complex and computationally intensive 
Kalman filter. Variational methods minimize the distance between the model and data over 
a time interval, subject to a dynamical or smoothness constraint and are also 
computationally intensive. This research directly benefits from the many studies done on 
earlier models which used a variety of techniques and evaluated the effectiveness of the 
assimilation. 

1. Sequential 

Nudging and Kalman filters are methods which attempt to move the model 
estimates of the flow fields towards observations. They differ, in that the Kalman filter can 
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Figure 2. Ratio of the non-linear terms to the geostrophic terms greater than 
10%.Blue is positive and red is negative, a) U components and b) V components. 
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estimate errors in the adjustment that is being made. Simply put, given x(t) as the model 
estimate of a model prognostic variable at time t and y as the observation at that point, then 
the model corrected estimate is: 

x(t)coirected = x(t) + a (y(t)-x(t)) (1.4), 

where a is a weight assigned to the observation. How a is defined and used is where the 
nudging method and the Kalman filter differ. Nudging, introduced into the meteorological 
field hyAnthes (1974), expands (1.4) into an equation written as: 

^=G + R(y-x) (1.5), 

where R is the nudging factor, x is the model estimate and y is the observation. G is the 
physics in the equation representing all the linear and non-linear terms for estimating the 
variable x. R can be written as a relaxation value, constant for all of space and time, such 
as is done for allowing for heat flux changes in a model, that is; R is equal to the inverse of 
several days. Complicating the equation slightly, R can be modified to include changes in 
space, accounting for known errors in the observation field. Additionally, a multiplicative 
term can be incorporated to account for the age of the observation. An observation that is 
old should not carry the same weight as a more recent observation of a neighboring point. 
As an example, Holland and Malanotte-Rizzoli (1989) use 

-(x2+y2)/L2R -Yt, 

a = ttQe ^ ^ (e ) (1.6) 

showing both a gaussian, (x^+y^)/L^^, based on the first Rossby radius of deformation 
(Lr), for spatial weighting, and an exponetial time decay, -yt, for temporal weighting. 

Different modifications of the basic equation (1.4) produces the Kalman filter 
equation: 
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^ =G + K(y-Hx) 



(1.7), 



where H is the matrix that denotes where there are observations and how the data has been 
interpolated onto a model grid and G, the model physics equal to Ax(t-l) +F, with A the 
transition matrix and F, the external forcing. Also, 

K = PH'’''[HPh'’''+R]‘^ (1-8), 

where P is the error covariance matnx for the model estimate equal to AP(t-l)A + Q (the 
process noise) and R is the error covariance matrix for the observations. H and R can be 
determined prior to a model run, but the matrix P can only be determined during the model 
run at specified time increments. Thus, the inversion of a portion of K is computational 
intensive to find the new estimate for x and makes the use of the Kalman filter impractical 
for a high resolution simulation. The filter does provide a way to continuously monitor the 
error in the solutions: Panalysis> for t, can be defined as equal to [ I - KH]P, where 
I is the identity matrix. 

There are two concerns when the filter is applied to non-linear problems. Although 
the filter is robust for linear problems, it is only approximate for non-linear ones. The 
second problem is the computationally intensive nature of the Kalman filter. Researchers 
have developed methods to overcome these problems, somewhat. One approach is the 
Extended Kalman filter (EKF) (Budgell, 1986), which approximates the error evaluation 
by piece-wise Unearization. The second simplification that has been done is to assume that 
P asymptotically approaches a steady state (Fukumori et al. 1993) and is calculated only 
once for equation 1.8. Fukumori and Malanotte-Rizzoli (1995) make another 
approximation by only assimilating the large-scale structure into the model. Thus, the 
Kalman gain (K) matrix is reduced in size and easier to handle. 
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2. Variational 



The adjoint method is the variational method most commonly used in oceanographic 
assimilation efforts. As an inverse problem formulated as an optimization problem, the 
method consists of minimizing a cost function (the model's numerical equations are the 
constraints) to find the minimum efficiently {Tziperman and Thacker^ 1989). A gradient 
based algorithm is used in which the gradient is computed using the adjoint equations of a 
model. Unknown model parameters are calculated by fitting model data to observations, 
with the addition of an optimal estimate of the observational data. Unlike the Kalman filter, 
this method is not directly related to the assimilation technique used in the present research 
and therefore, the mathematics of the method will be omitted. The method works as 
follows. Using some guess for the initial values, estimates of the prognostic variables at 
t+1 are computed and used in the formation a cost function, J, measuring how well the 
estimates fit the observed values. In an inverse and interactive fashion, the cost function, 
J, is then used to adjust the initial values such that J is minimized. Obviously, running the 
model forward, then backward through the use of its adjoint, at least doubles the cost of the 
computation, if not more. Thus, the procedure is undesirable for a high resolution global 
ocean model with over 30 million values to be calculated. Sirkes et al. (1996) used an 
adjoint technique applied to a coarse grid global ocean model in an attempt to improve 
estimates of eddy parameterizations and to produce a better initial ocean model for coupled 
ocean-atmosphere climate simulations. Although the results are encouraging, there are still 
difficulties in the formulation relating to conditioning of the dynamic constraints. 
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D. OVERVIEW OF ALTIMETRY: OBSERVATIONAL DATA 



FOR ASSIMILATION 

The ERS-1/2 {Le Traon et al. 1995) and T/P (Fu et al. 1994) satellites are now 
providing researchers with high quality and long altimetric time series for use as the 
observational field to adjust OCGM prognostic fields. In this research, data from both 
these instruments has been combined to form a joint sea surface height field which is to be 
assimilated. 

An altimeter measures the distance (ha) between the surface of the earth and the 
location of the satellite. This is done by having the satellite emit a pulse of radiation and 
then computing ha by using the time it takes the pulse reflected from the earth's surface to 
return to the satellite. The height of the sea level (SSH) is determined by subtracting ha 
from the total height the satellite is above an orbital reference ellipsoid (H) minus the height 
of the earth's geoid (hg) (also referred to as a level of no motion). Additionally, corrections 
(hatm) are included in the SSH estimate to account for atmospheric conditions (both dry and 
wet tropospheric); it takes longer for radiation to travel during wet conditions than dry, in 
addition to the effect the atmospheric pressure (hinv) has on the sea surface (known as the 
inverse barometer effect). If the researcher is interested in the general circulation of the 
ocean, then any contribution the tides (htides) make on SSH needs to be removed also. 
Thus, the total equation for determining SSH from what the satellite measures is: 



SSH (i,t) = H(i,t) - ha (i,t) - hg(i) - hatm(i4) ■ hinv(i»t) ■ htides(i4) (1-9), 



where i is a spatial location and t is the time of the measurement. 

Unfortunately, for wavelengths shorter than about 2500 km {Rapp et al 1996), hg 

is not known accurately and complicates how the data may be used. Equation 1.9 can be 
rewritten as follows with SSH separated into a mean ( SSH ) and a residual SSH' and the 
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two static terms, SSH and hg combined; 



SSH'(i,t) + [hg(i)+ SSH (i)]— H(i,t) - ha(i,t) - hatm(t»0 * hinvOiO' bndesCbO (1.10). 

As the equation is shown in 1.10, the terms on the right are all known at the time the 
satellite makes its measurement. H is determined from orbital mechanics using a set of 
either lasers or other tracking devices located on a satellite to track its position (for 
Topex/Poseidon, it is the DORIS, Dual-Doppler Tracking System Receiver, which receives 
signals from ground stations for satellite tracking, for ERS-1/2 it is PRARE, Precise Range 
and Range-rate Equipment^). The atmospheric corrections (hatm) are determined either 
from a model and/or the use of the dual frequency altimeter on T/P, and the microwave 
radiometer on T/P and ERS-1/2. The tidal correction (htides) for both altimeters is 
determined from an empirical model (CSR 3.0, Eanes and Bettadpur, 1995). And last, ha 
is the measurement made by the satellite (averaged to 1 Hz from 50 Hz). Thus all the terms 
are accounted for except for the two static terms. The fact that these satellites fly in orbits 
which repeatedly sample the same ground track, can be used to determine the value of the 
combined static terms by finding the mean over time of [H - ha - hatm - hinv - htides] at 
specific spatial locations. Because there are two static terms in the equation 1.10 and 
neither are known to the accuracy needed for use in studying the ocean, only the mean of 
the two terms together can be found and removed, leaving only the sea level residual or 
anomaly, SSH'. Specifics on the data processing required to use altimetric data can be 
found throughout the literature on satellite altimetry and ocean circulation (e.g. Tokmakian 
and Challenor, 1993). 

One additional wrinkle must be taken care of before the altimeter data can be used 
for assimilation. The assimilation procedure needs to use SSH and not the residual. Thus 

^ The PRARE failed on ERS-1 and the orbit was mostly determined from the use of Laser tracking of the 
satellite from the ground. 
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a mean sea level field must be added to the anomaly measurement prior to being 
incorporated into the model. For the purposes of this research, a mean sea level field was 
used from a prior run of the model. Analysis by Rapp et al. (1996) shows that the 1/4° 
POCM model mean field closely proximates the mean sea level as estimated empirically 
from observational data (mainly T/P data) at long wave lengths. 
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II. ASSIMILATION TECHNIQUE 



The emphasis of this dissertation is to assimilate altimeter SSH data into a high 
resolution PE model and due to the computation costs of the other methods, the 
assimilation options are restricted to that of nudging. Two experiments have been 
conducted. The first one is referred to as a twin experiment. For this experiment, two 
model runs are made. The first is the control run and the second is the twin, a simulation 
which uses simulated SSH data from the control run in an attempt to reproduce the 
conditions seen in the control run. The second experiment uses real altimeter data to force 
the model more towards what observations show. Unfortunately, the nudging technique 
does not include a way to assess the error in the estimated field produced by the 
assimilation as it is compared to the observational field of data; thus some simple statistical 
tests are also described that are used to evaluate the success of the experiments by 
comparing the results to the control runs or to the observational data. 

A. NUDGING SPECIFICS 

The specific equations for the assimilation follow. The height field is assimilated in 
a straight forward manner, namely 

^ = G + R*(T|0-Tim) (2.1), 

where m denotes model data and o denotes observational data, which in this case is the 
altimeter anomaly height plus a mean height of a previous model run. G is the model 
physics for this particular equation and R is the relaxation coefficient. R is a function of the 
error in the observation field times a temporal factor times a spatial factor (Anthes, 1974). 
For example, R at grid point ij is 
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( 2 . 2 ), 



-5t/T -(x2+y2)/LR2 
Rjj = £ e e 



where e is an error estimate of the model data point. In this case, e is related to 



(little variance), the point will be nudged weakly. When it is large, the model SSH will be 
pushed harder towards the observed value. This term includes observational statistics 
when nudging the model fields and is somewhat like a Kalman filter which continually 
updates its error covariance matrix (unlike this one, which is static). The observational 
SSH provides the two-dimensional set of predetermined e values, calculated prior to the 
assimilation run. The last two terms weight the observation with respect to time and space, 
an exponential decay in time and a gaussian distribution in space on the order of the size of 
the Rossby radius of deformation (Lr). Daily, the SSH altimetry (either synthetic or 
observed) data set is gridded to the model grid and a decay factor, 5t, is aged. For 
example, a new observation has a decay weight of 1 (5t = 0) and an observation 5 days old 
has a weight of 0.37, if T is equal also to 5 days. Therefore, as an observation ages over 
time, it is given less weight. 

There are several ways to relate to surface height signature to subsurface quantities. 
Hurlburt (1986) and Carnes, et al. (1990) use empirical orthogonal functions (EOFs) of 
historic temperature profiles to relate the surface height and temperature to subsurface 
fields. In the EOF method, the vertical profile of temperature is decomposed into 
orthogonal modes allowing for temperature estimates to be made based on its variance, 
giving as the EOF equation: 



<5h(i,j)>2 



the normalized variance of SSH at a point. Thus, when the value is small 





(2.3), 



n 
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where 0 is the new estimate of temperature at z, <T> is the mean temperature for that 
point, the eigenvalues, and fn(z) is the value of the eigenvector at z. To relate SSH 
to these EOFs, Carnes, et al. fit values of measured SSH to the eigenvalues over all space, 
such that X.n = bo + biT) + b2T|2. With this new estimate of the eigenvalues, the variance 
matrix is more like the covariance of <T*t|> rather than <T*T>. The vertical structure f(z) 
is assumed to be constant over the whole region for monthly periods. SSH is a measure of 
the change in the barotropic component of the flow and thus, the corrections made during 
assimilation runs are just scaling factors relating to the change in SSH, and are barotropic 
in nature. But in reality, especially on a global scale, adjustments need to be made also to 
salinity. 

Another approach has been taken by Haines (1991) and applied to a PE model by 
Cooper and Haines (1996). This is a dynamic method which adjusts the density field up 
or down while conserving potential vorticity. This method requires iteration at each grid 
point, making adjustments, in theory, to both the temperature and salinity fields, until the 
correct surface pressure is achieved and potential vorticity is conserved. A side note. 
Cooper and Haines use a simplified equation of state: p = 1029. + 0.25T, which makes it 
necessary only to make changes in temperature, rather than both temperature and salinity. 
The method can be computationally intensive and thus not appropriate for a high resolution 
global model. 

Mellor and Ezer (1991) found, using a primitive equation model, that predetermined 
correlations between sea surface height anomalies and subsurface temperatures and 
salinities were required to translate surface information to the subsurface layers. This is the 
method chosen for this research. Here, correlations are found between SSH and T for each 
layer at each grid point and between SSH and S for each layer at each grid point. This type 
of correlation is more of a baroclinic correction, although, if all the layers have the same 
correlation between SSH and T and S, then it becomes a barotropic correction. These 
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relationships are used to create an estimate of the observed subsurface temperatures (or 
salinities) which correct the simulated fields towards observations. The basic equation for 
including the correction term in the temperature equation is 



-^ = G + R* FT(TlO-Tin’) 



(2.4), 



where is an additional term correlating SSH to the subsurface temperature. By replacing 
T with S a corresponding equation for salinity is defined. The correlation functionals for 
each model grid point are written as: 



Few global temperature and salinity data sets exist that include a temporal 
component. Therefore, the twin experiment (see Chapter III) uses statistics from the 
control run to determine these functionals and the Levitus 1994 (Levitus and Boyer, 1994) 
monthly climatology of T and S has been used in the experiment using observational 
altimetry (see Chapter IV). The covariance matrix for the experiment using real altimeter 
observations is computed using the Levitus monthly temperature at a grid point with a 
monthly SSH field computed from the altimeter data set. This produces a three 
dimensional set of functionals that is used in the nudging of subsurface fields. It might 
also be possible to create these global functionals using data from different climatology sets 
for different oceans, although it would be necessary to examine the compatibility of one 
data set to another. For example, the Optimum Thermal Interpolation System (OTIS) 
temperature data set, produced by the US Navy in the North Pacific and North Atlantic 
oceans, might be used instead of the Levitus 1994 monthly climatology. The method can 
also be modified to assimilate actual temperature or salinity values. 



<ST5ti> 

<5t|2> 



and F^= 



<5S5r|> 

<6t|2> 



(2.5, 2.6). 



18 



B. DESCRIPTION OF STATISTICAL TESTS USED FOR 



EVALUATION 

Basic statistical methods can provide useful ways to evaluate large data sets. Three 
tests are employed to describe the major differences or simularities between model runs, the 
Student-t, chi-square, and F tests. With a set of paired runs (control and assimilation), 
statistical tests such as these are very appropriate in determining the effectiveness of the 
propagation of altimeter data into the deeper levels of the model. 

1. Student-T Test 



The Student-t test is used to determine if there is a significant difference in the mean 
of the distribution of one prognostic variable in a given region for the two model runs as 
shown in the following explanation. Let a/ be a data value at a location from the control 
run and hi be a data value for the same variable and location from the assimilation mn. Dj 
represents the difference in the two values (D/ = a/ - b/ ). It is assumed that the 
differences (Dj) are distributed about a mean ^ and that the deviations are normally and 
independently distributed about this mean. The Student-t distribution tests the null 
hypothesis that jj=^. The quantity. 



(<D> - p.) 
Sd 



( 2 . 7 ), 



follows the Student-t distribution with (n-1) degrees of freedom, where <D> is the mean 
over all pairs of the D/ , sd is the estimate of the standard deviation of the differences: 




(Pi - <D>)2 
V n(n-l) 



( 2 . 8 ), 



19 



and n = number of pairs (Snedecor and Cochran, 1989). For this paper, n is always 
greater than 120 and the probability value from the table for (n-1) = co can be used. That 
is, if t > 3.29, there is less than a 0.1% chance that the null hypothesis (|X=0 ) is true. In 
other words, the two populations have different means. 

2.Chi-Square Test , 



The goodness of fit test is used to examine the similarity of the distribution of 
the data from the two sets of data. Let o/ = the number of observed values (or values from 
the assimilation model run) that fall into bin i and ef is the number of expected values (or 
values from the control run) that fall into bin i. Then, 



= S 



(o/ - e/ )2 

ej 



(2.9), 



where the sum is over all bins. The null hypothesis is: the observations are from the same 
population as the expected population. For the tests applied to the data presented in this 
paper, the number of bins are in the neighborhood of 25, some required more and some 
less. In general, a conservative approach is to determine if the distributions are similar or 
not. For example, if the number of bins (n) equal 25 (degrees of freedom: df = n-1) and 
the calculated greater than 47 (determined from tables, see Snedecor and 

Cochran, 1989), then there is a less than 0.5 % chance that the two distributions are the 
same. The advantage of the test over the Student-t test, is that the the yp- does not 
expect the data to be distributed normally and to be independent. 
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3. The F Test 



Related to the y} test is the F test. The F test, when used to compare two 
populations, determines if the null hypothesis is true: Si^ and S 2 ^, sample variances from 
independent samples, have true variance Given that F = then the null 

hypothesis will be satisfied if the observed F < Fo.o25- Fo.025 is determined from tables 
such as is found in Snedecor and Cochran, 1989. For large degrees of freedom, F 
approaches 1 when the null hypothesis is true. 

Theoretically, these statistical tests should show that fields from a twin experiment 
should be similar and that when real altimeter data is assimilated, the fields of the control 
run and the assimilation run should be statistically different. 
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III. TWIN EXPERIMENT 



A twin experiment was conducted to determine if the assimilation process could 
reproduce a known mode state (i.e. the control state) when initialized from a different 
starting location. The twin experiment, thus, is a test that has known initial and final states 
which can be used to evaluate the effectiveness of an assimilation process. To simplify the 
text when describing the results, the control run will be referred to as run "C", the 
assimilation run (produced using synthetic altimeter data from the control run) as run 'T" 
and the run which starts at the same point as "T" but without assimilation is run "I" for 
initial run (or field). The rest of the chapter gives the details of the test and the results. A 
success for the twin experiment is if, at the end of the T run, the prognostic and diagnostic 
(such as potential vorticity) fields resemble the control run, C. Analysis will not be done in 
this section which shows if the model, itself, is simulating a realistic ocean. 

A. SETUP 

The control model run is from the simulation referred to as the Semtner/Chervin or 
POCM_4B. The year 1992 from this simulation will be the control for the twin experiment. 
An evaluation of this model run is described in detail in Stammer et al. (1996). The model 
was forced with wind stress fields created from nine years (1987-1995) of 10 meter 
ECMWF wind velocity fields and the Barnier et al. (1994) monthly climatological heat 
flux. The twin experiment model run (T) is in the exact same configuration as the control 
run (C), with the same forcing, with the exception of the inclusion of numerical code to 
handle the assimilation of simulated altimeter data and an associated vertical structure. The 
T (or assimilation) run is initialized with a three dimensional field for day 362 of 1987, a 
statistically different field than that which is used to initialize the C run (day 365.9 of 
1991). The twin experiment is to determine if, by starting at a different initial field, the 
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model can be nudged towards simulating the ocean's circulation as seen in the control run. 
The twin experiment was conducted for a period of one year. 

Figure 3a shows the initial SSH field for run T and while 3b shows a similar SSH 
field for run C. The difference between the two fields is shown in Figure 3c. Large 
differences are clear in the North Pacific, where the SSH field from the initial day (Figure 
3a) is much higher than from the initial field from the C run. In the Antarctic Circumpolar 
Current (ACC) and in the North Atlantic, the differences appear more on the mesoscale, 
rather than on a basin wide scale. 

1. Synthetic Altimeter Observations and Associated Vertical 

Functionals 

A simulated sampling pattern associated with the 35 day repeat period of ERS-1 in 
both time and space was created and then used to sample the SSH (x\) of run C for the year 
1992, every three days, creating a synthetic SSH assimilation field for every day of 1992. 
The ERS-1 sampling is replicated in both time and space for this experiment. Figure 4 
shows an area from the North Atlantic which shows the grid points that are sampled in a 
ten day period. At 60°N, the track separation is at 39km, nearly the same as the grid 
spacing of the OCGM. Near the equator, the altimeter sampling is about 50-60km, about 
the sampling of every other model grid point. These fields are the observations to be 
incorporated into the model. Figure 5 shows the corresponding synthetic observations for 
early in December of 1992. It shows the dense coverage of the satellite track coverage of 
the model grid. 

Next, equations 2.5 and 2.6 are used to form the vertical functionals which 
associated SSH with subsurface temperature and salinity. The covariance quantities: 
<5T5 ti>, <5S5i)>, and <5ii2> are formed (for this experiment) using three years of 
previous model output covering the years 1987 through 1989. Because all the observations 
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a) Initial Field T Run 
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b) Initial Field C Run 
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Figure 3. a) Initial SSH field for T run (December 28, 1987). Units in centimeters, b) 
Initial SSH field for C run (December 28, 1992). c) Difference between T field and C field. 
Scale in units of centimeters. 
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and associated vertical structures are from runs of the same model, their dynamical 
properties are the same and the Helds are dynamically consistent. Averages of the vertical 

functionals over 20' squares have been formed and normalized in each layer by 

F tempk ^ standard deviation of layer 1 , , . . , 

: *Tr — = — r * 1 — j — j"j — ft — » where k is the layer, to give a 

maximum value in k standard deviation of k j » & 

relative picture of how the vertical structure in either temperature (Figure 6a) or salinity 

(Figure 6b) are nudged. For example, in Figure 6a, in the 20' box centered at 370'E, 5°S, 

temperature is nudged strongly at the surface and very little at depth, while at 250'E, 45°S, 

there is little relationship between the vertical temperature variations and the height 

variations. 




Figure 4. Example of the sampling grid over ten days of the model grid for creating 



synthetic altimetric SSH fields for assimilation during the T run. 
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Figure 5. SSH field corresponding to the sampling field in Figure 4. Units in centimeters. 

2. Additional CoefHcients in Nudging Expression 

Included in the nudging equation (Equations 2. 1 and 2.2) are two additional terms 
which need to be calculated to contribute the weighting term associated with the nudging 
scheme. The first term is a weight, e, associated with the normalized variance of SSH at a 
grid point with the underlying relaxation coefficient at three days. This term is created by 
using the variance of SSH over three years (1987-1989). Also included in the weight is an 
exponential which extends the data, sampled at a grid point, over a wider area. This term is 
a gaussian term with the Rossby radius of deformation (Lr) as the controlling term. Lr is 
determined also from the model data using a 3 year mean. 
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Figure 6. a) 20°x20° mean, normalized functional for temperature produced from equation 
2.5 in Chapter II. Y axis of each square refers to the model layer, with 1 being at the 
surface. Blank squares are over land. 
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B. STATISTICAL COMPARISONS 



1. RMS Differences 

Root mean square (RMS) differences of SSH, temperature, and salinity over the 
global oceans and individual basins over time is one measure of how a model is adjusting 
the assimilation process. Figure 7 shows the time series of RMS differences between the 
SSH of the T and C runs over the first 250 days of 1992. The differences are normalized 
by the maximum difference of each time series. This allows for the comparison of how 
each region is adjusting, in addition to the global adjustment. The global difference has 
dropped to about 50% of the initial value between the two runs, with the South Atlantic 
basin dropping the quickest. Much more significant adjustments are occurring where the 
largest differences are seen in Figure 3c, in the Indian and the North Pacific, which have 
reduced the error to the order of 25% and 30%. These values are similar to both Oschilles 
and Willibrand (1996) and Cooper and Haines (1996) which use only regional high 
resolution PE models, with different assimilation schemes and somewhat idealized forcing 
(climatology instead of daily varying wind stress and monthly varying heat flux). These 
curves show that the present method is having some success, especially considering that 
the model is not being adjusted over large areas of the world's ocean where the variability 
is changing little or not at all (related to the e term of the nudging coefficient). 

Similarly, the time series of the RMS difference in temperature and salinity is 
shown in Figure 8a and b. The plots show the global RMS differences in temperature and 
salinity over time for the top twelve layers of the model, down to approximately 1200 
meters normalized by the maximum difference seen in each layer. As expected, the largest 
change is seen in the near surface layers. But even in layer 12, the temperature difference 
has been reduced by 30% by day 250. The decrease and then the increase (around day 
100) in the temperature difference in layer one is due to the adjustment the model is making 
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to the heat flux term, rather than to the assimilation process. This can be compared to the 
salinity RMS differences (Figure 8b) where the difference drops and then stays low. 




Figure 7. Plots of time series of RMS differences between T run fields and C run fields of 
SSH. Solid black line represents the global RMS differences over time, solid red 
correspond to the North Atlantic, dotted red line: South Atlantic, solid blue line: North 
Pacific, dotted blue line: South Pacific, and green line to the Indian. Values have been 
normalized by the maximum RMS difference calculated for the specific basin or globe. 
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rms diff temp 1-12 layers (0 — > 1 160meters) 




Figure 8. Same as for figure 7, but showing the RMS difference in a) temperature and b) 
salinity for each layer. Values have been normalized by maximum RMS value seen in each 
layer. The numbers on plot identify the layer. 
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2. Statistical Tests 



The two statistical tests, the student-T test and the F-test, are used for a quick 
evaluation of the distributions of various model prognostic variables for different regions. 
It is important to note the assimilation of data using the nudging method will not 
immediately adjust the values of the model's variables to the observational values. Rather, 
the model is gently pushed towards what is observed. Changes in the deeper layers will 
also be difficult to see immediately, mainly because over large parts of the ocean, the deep 
ocean is decoupled from the surface. Thus, these simple tests are useful to determine if, in 
time, the monthly averages are approaching the values of the control field. 

Figure 9a-h show two sets of student-t values for temperature from the twin 
experiment. Using equation 2.7, the t-test value is calculated between the various mean 
fields for the globe and for various regions: North and South Atlantic, North and South 
Pacific, Indian, tropical oceans and the Southern Ocean. Some of these regions overlap. 
The vertical scale is the t-test value with any value above 1.64, statistically different. Along 
the horizontal axis, the t-test value for SSH is at x=l and for the 20 layers of temperature, 
the values are at x = 2-21. On the left side (Figures 9a-l,9b-l,9c-l, etc.), the plots show t- 
values computed between different months of the control run, December 1987 and January 
1992 (thin dash-dot line), and between December 1987 and December 1992 (bold dash-dot 
line) to establish a baseline of how the model differs between years without any 
assimilation. Note the values are very similar between the two lines, especially below layer 
three. Statistically, these plots suggest there is no difference in the means of temperature 
between the initial field when it is compared to January of 1992 and when compared to 
December of 1992. These plots also suggest what differences to expect between two years 
of the model when no assimilation is performed. 
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Figure 9. Plots showing the t-test values calculated between different model fields for the 
globe and for different basins. A pair of panels comprise each figure. Y-axis is the t-test 
value, 1.64 is the significance level for a null hypothesis. X-axis is read as follows: x = 1, 
t-test value for SSH, x=2-21, t-test value for temperature, levels 1-20. There are six areas: 
a) global values b) North Atlantic, c) South Atlantic, d) North Pacific, e) South Pacific, f) 
Indian, g) Tropics, h) Southern Ocean. There are two panels per figure letter. The left 
panel represents the t-test values between a field from December 1987 (I field, initial field 
of T run) and January 1992-C run (thin dashed line) and December 1987 and December 
1992 (bold dashed line). On the right, the t-test values have been computed between the 
Januarys of the T and C runs (thin line) and the Decembers of the T and C runs (bold line). 
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The right set of panels show corresponding plots of t-values between the T and C 
runs (Figures 9a-2,9b-2,9c-2, etc.) for January 1992 (thin line, early in the simulation) and 
for December 1992 (bold line; at the end of run). Since the T run has been initialized with 
the December 87 field, it is expected that the thin lines on Figure 9a2, b2, etc. are close to 
the thin lines on Figure 9al, bl, etc. If the two plots are overlayed, one on top of the 
other, they are quite similar. This means that in the early part of the assimilation process, 
the mean of the temperature field has changed very little. When comparing the two lines in 
the right panels (Figure 9a2, b2, etc.), it is seen that in all regions, except the tropics 
(Figure 9gl and 2), there is a reduction in the t-test values from the beginning of the year 
(thin line) to the end of the year (bold solid line), much greater than seen in the t-test value 
between the years in the control run (Figures 9a 1, bl, etc.). These sets of plots also 
suggest the differences in the means near the surface are mostly due to differences in the 
wind stress fields forcing the model, exhibited by the differences in the t-values shown in 
the panels on the left side because only the wind stress is significantly different in the two 
runs. 

The results shown in Figure 9 are consistent with the RMS differences that shown 
in Figure 8, that is, the assimilation is affecting the temperatures down to about 1000 
meters, with very little change below. What these plots show is that there has been greater 
change in the model's temperature and SSH fields when assimilation is performed, then 
when it is not; and that change is towards the control field, rather than away because the 
values are getting smaller, consistently. 

Figure lOa-h show a similar set of t-test values for the salinity field. Similarly to 
the patterns for temperature, the pattern of reduced t-test values (Figure 10a2, b2, etc. ) is 
seen for layers 1 through 10 for the t-test over those values computed between two years of 
the control simulation (Figure lOal, bl,etc.). There seems to be little change in the pattern 
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Figure 10. Same as 9, but representing t-test values computed from salinity fields. X=l- 
>20 corresponds to model levels l->20. 



in the North Atlantic (Figure 10b) because aU the calculations show that difference between 
the December 1992 T and C runs (bold line in Figure 10b2) is very similar to the t-test 
values calculated for the comparison between the December 1987 fields and December 
1992 fields of the C run (dashed bold line in Figure lObl) when the two years are very 
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different. The Southern Ocean (Figure lOh) exhibits a much different picture. Here it 
appears that the model is making significant modifications to the salinity field to adjust it 
towards the control field. 

As stated in Chapter II B3, the F-test is a technique to examine the differences in 
variance between two sets of data. Figures 1 la-h show the computed F values for the total 
globe and various regions for SSH and temperature in a similar fashion to that done for the 
t-test. The vertical scale is the F-test value with the significance level at 1 .0. The left 
panels of the figures are the values computed between different years and months of the 
control field (C Dec. 1987 - C Jan 1992 and C Dec 1987 - C Dec. 1992, dash-dot lines); 
and the right panels are the C-T comparisons (T Jan 1992 - C Jan 1992 and T Dec 1992 
and C Dec 1992). Again, as in the Student-t calculations, the comparison of control fields 
do not change much between the different months, while between mns T and C (solid 
lines), the F-test computed values have changed and in most cases, are adjusting towards 
1. Figure 12a-h show the same thing, but for salinity. It is much harder to detect any 
differences in the variability in the temperature and salinity patterns from these figures, 
probably only in the top layers of the model. But, this should not be that surprising, given 
that the observations are from the control run and one would expect similar variability for 
both simulations. 
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Figure 11. Same as 9, except the Y-axis represents the results of F-test computations, with 
the significance level at 1.0. 
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Figure 12. Same as 1 1, except representing F-test values computed from salinity fields. X 
=l->20 corresponds to model levels l->20. 
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C. POTENTIAL VORTICITY 



f 3p 

A simple, linear definition of potential vorticity (PV), Q = — where p© is the 

Po 

average density, f is the Coriolis parameter and z is the depth, can be used to examine if 
changes have occurred in the density structure of the T run so it is simulating the control 
run rather than the I run. Q has been calculated along several density surfaces, Oq = 26.85, 
27.0, and 27.3 with a thickness of 0.15 kgm'^ (corresponding to how model layers are 
defined in an isopycnal model. New et al. (1995) for September 5, 1992 for both the C and 
T runs and for December 28, 1987 for run C, the initial field for run T. Table 1 
summarizes the results by listing the RMS differences between the two September fields 
and between each of the September fields and the initial field for the assimilation run in the 
left three columns. The right three columns are comparing similar RMS differences but 
using fields from the early part of the experiment (January 2, 1992). RMS differences for 
both potential vorticity and thickness of the a layer are listed in the table. The table 
consistently shows, using the RMS differences, that late in the experiment (September 
1992), the assimilation run (T) is reproducing the potential vorticity fields calculated from 
the control (C) run at all three a levels (lines 1 ,4, and 7, column 4) much better than early 
in the run (lines 1,4, and 7, column 6). The first 3 columns also show that the RMS 
difference values produced when comparing run T and run C to the initial field are 
approximately the same for both potential vorticity and layer thickness. In the right most 
three columns, it is the RMS differences between the field from a early point in the 
assimilation compared to the initial field which correspond the best. This is as expected; 
the early period of run T resembles the initial field, then slowly is nudged more closely to 
the control run. 
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FIELDS 

End of twin expr. 


RMS 

diff 

OelO 


RMS 
Az 
. (m) 


FIELDS 

Begin, of twin expr. 


RMS 

diff 

QelO 


RMS 

Az 

(m) 




Oq = 26*85 






Oq “ 26.85 






1 


Sep 5 92 (C) - Sep 5 (T) 92 


2.41 


35.51 


Jan 2 92(C) - Jan 02 (T) 92 


3.92 


55.08 


2 


Dec 29 87 (I) - Sep 5 (T) 92 


165.22 


25.97 


Dec 29 87 (I) - Jan 2 (T) 92 


1.76 


22.21 


3 


Dec 29 87 d) - Sep 5 (C) 92 


165.24 


35.22 


Dec 29 87 (D - Jan 2 (C) 92 


4.07 


55.62 




Oo = 27.00 






Oo = 27.00 






4 


Sep 5 92 (C) - Sep 5 (T) 92 


1.75 


37.99 


Jan 2 92(C) - Jan 2 (T) 92 


3.35 


52.94 


5 


Dec 29 87 d) - Sep 5 (T) 92 


113.15 


29.73 


Dec 29 87 a) - Jan 2 (T) 92 


1.85 


22.00 


6 


Dec 29 87 (I) - Sep 5 (C) 92 


113.16 


40.76 


Dec 29 87 (I) - Jan 2 (C) 92 


3.57 


53.84 




Oo = 27.30 






Oo = 27.30 






7 


Sep 5 (C) - Sep 5 (T) 92 


1.69 


33.63 


Jan 2 (C) - Jan 2 (T) 92 


3.11 


43.26 


8 


Dec 29 87 Q - Sep 5 (T) 92 


2.89 


30.00 


Dec 29 87 a)-Jan2(T)92 


2.22 


19.87 


9 


Dec 29 87 H) - Sep 5 (C) 92 


3.16 


37.35 


Dec 29 87 (I) -Jan 2(C) 92 


3.00 


43.42 



Table 1. RMS differences in Q and layer depth (Az) for various pairs of days of the model 
runs. September 5, 1992 is near the end of the twin experiment. January 2, 1992 is near 
the beginning of the experiment. December 29, 1987 is the field used to initialize the 
assimilation run of the twin experiment. 



Figures 13, 14, and 15 a-c give plan views of the potential vorticity on the ao=27.0 
isopycnal while d and e show the differences between the September 5, 1992 C and T 
fields and the September fields from the T mn and the initial field near Dec. 29, 1995. The 
changes in PV for the different runs (C and T) over the North Atlantic (Figure 13) are most 
obvious in the very north of the model domain. The 27.0 isopycnal surface in the T run 
(Figure 13a) no longer outcrops along a diagonal from 45°N, 300°E to 55"N, 350°E as 
seen in the initial field (Figure 13c), but now is shown to have a similar PV value as for the 
C run (Figure 13b). Similarly, the lowest PV values along this isopycnal are similar to the 
C run, rather than the initial run. Figure 14 shows a similar view of the North Pacific. The 
band of high PV starting at about 40°N,160°E and continuing to about 190°E and spreading 
southwards is more representative of the C run (Figure 14b) than of the initial field (Figure 



41 



14c). In the tropical Pacific (Figure 15), the PV remains similar between the three runs, 
except in the eastern Pacific, where blobs of higher PV water are found in the T run similar 
to those seen in the C run. 



Twin Sept. 5, 1992 Control Sept. 5, 1992 




Initial Dec. 29, 1987 




Twin - Control Sept. 5, 1992 Twin Sept. 5, 1992 - InitiaI(Dec. 29, 1987) 




280 300 320 340 360 380 280 300 320 340 360 380 



Figure 13. PV on isopycnal surface = 27.0 in the North Atlantic, a) PV calculated from T 
run, day Sept.5, 1992. b) same but for C run, c) Initial field from T run. d) difference 
between T and C PV fields, e) difference between T and I PV fields. Units are <Qel0>. 
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Twin Sept. 5, 1992 



Control Sept. 5, 1992 




Initial Dec. 29, 1987 




Twin - Control Sept. 5, 1992 Twin Sept. 5, 1992 - Initial(Dec. 29, 1987) 

1.5 










200 



250 




Figure 14. Same as Figure 13, except for North Pacific. 
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Twin Sept. 5, 1992 Control Sept, 5, 1992 




Twin - Control Sept. 5, 1992 Twin Sept. 5, 1992 - Initial(Dee. 29, 1987) 




Figure 15. Same as Figure 13, expect for tropical Pacific. 
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D. REGIONAL ANALYSIS 



1. Southern Ocean 



Three meridional sections across the Southern Ocean have been examined. From 
east to west, the sections taken are across the full width of the Southern Ocean at 20°E 
(South Africa), 150°E (Australia), and at 290°E (Drake Passage). The meridional sections 
extend from the model's southern point at 75°S to 35°S. Volume transports across the 
three sections are listed in Table 2, calculated from instaneous fields from Sept. 1992 and 
the initial field from 1988. 



Section 


Initial 


Control 


Twin 


S. Africa 20°E 


162 SV 


160 SV 


167 SV 


Australia 150°E 


164 SV 


166 SV 


172 SV 


Drake 290°E 


166 SV 


161 SV 


167 SV 



Table 2. Volume transports across three zonal sections of the Southern 
Ocean at 20°, 150°, and 290°E from 75°S to 35°S. Units are in Sverdrups, 
1SV= 10e6 mh-K 



The mean of all the transports is 165 +/- 3.5 SV, with the highest transport values 
found in the twin experiment across the Australia section (172 SV). After the assimilation, 
the run T has produced transports higher than run C by around 6 to 7 SV for all three 
sections. Both have consistent values between the South Africa section and the Drake 
Passage. Figure 16 shows the SSH from snapshots of the model on Sept. 5, 1992 for the 
runs T and C and the initial run I (Dec. 29, 1987). For the Drake Passage (Figure 16c) and 
for the meridian south of Australia (Figure 16b) the adjustment of SSH in run T (dashed 
line) towards the control run (solid line) can be seen. The red line represents the relative 
strength of the nudging coefficient at the surface. All the differences between the initial 
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field (dotted lines) are not due to the assimilation of SSH. Some of the differences are due 
to the application of a different wind stress field. However, when changes appear 
consistently over a group of latitudes such as in Figure 16b between 65°S and 55°S and 
around 43°S, the change is mostly due to a large change in the underlying vertical structure. 



SSH 20E 




latitude 



Figure 16. Meridional plots (black lines) of SSH in centimeters at 20°E, 150°E, and 290°E 
between 75°S and 35°S. Solid line is from C run, dashed line is from T run, and dotted 
line is the initial state from I. Red line is the relative weight associated with the nudging 



term. 



Figure 17, representing the differences in density between the snapshots of run T and the 
initial field (I) of the three meridional sections, shows that at 150°E (Figure 17b) the water 
has become less dense between about layers 3 and 12 (60 meters to 1000 meters, 
approximately) between 65° and 55°S, leading to an increase in the height field in the 
assimilation run (dashed line in Figure 16b), similar to the SSH of the control run (solid 
line, in Figure 16b). In contrast, between 43°S and 35°S, the density has increased with 
the assimilation process and thus, the SSH signal is reduced. 




Figure 17. Plots of differences in aobetween Sept. 5, 1992 T field and Initial field along 
the meridians at a) 20°E, b) 150°E, and c) 290°E between 75°S and 35°S. Scale in units of 
kgm‘3. X axis is model layer. 
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At 20°E (Figure 16a) , the large SSH signal produced in the Aguhlas retroflection 
region at -40°S, is not being adjusted as expected by the assimilation. Figure 18 shows the 
density differences between the T and C fields for Sept. 5, 1992. For the meridional 
section in density at 20°S (Figure 18a), the T run is substantially denser than run C (and 
thus, the SSH signal will be lower). Why the assimilation is not adjusting its deeper layers 
to the correct density can be found by examining the functionals (see Eq.2.5 and 2.6) 
which relate the surface heights to subsurface values. Figure 19 (temperature) and Figure 
20 (salinity) show the functionals for the three meridians across the Southern Ocean. At 




Figure 18. Same as Figure 17, except for differences in <Jo between Sept. 5, 1992 T and C 
fields. 
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approximately 40°S in Figure 19a, the functional for both temperature and salinity is much 
smaller than the values for the function at the same location in Figure 19b (for 150°W, 
south of Australia). This would suggest that it is the adjustment of the deep layers that is 
making the most change in the sea surface signal at these locations. 



Temperature functional 20E 




Temperature functional 150E 




Temperature functional 290E 




Figure 19. Meridional plots of temperature functional at 20°E, 150°E, and 290°E between 
75°S and 35°S. Y axis is model layer. Scale is in units of Cm’*. 
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Salinity functional 20E 




Salinity functional 150E 





Figure 20. Same as 19, but for salinity functional. Scale is in units of S m-^ (S is also 
defined as a practical salinity scale, PSS). 



2. Tropical Ocean 

A meridional section along 156°E from 5°S to 5°N (Figure 21a-j) shows the U 
velocity component for three fields. Sept 1992 T run (21a) and C (21b) run and the initial 
field (21c). Figures 21d-f , g-i are the corresponding temperature and salinity fields. The 
assimilation has adjusted the temperature field only slighdy, while the salinity field is 
much fresher in both the T and C simulations at a depth of about 80 to 200 meters (layers 
4-7). This freshening is reflected at the sea sinface, near -4°S and 2°N, (Figure 21j) with 
an increase in the SSH value in the T and C fields (dashed and solid lines) over the value 
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from the initial field, I (dotted line). The red line on Figure 21 j shows the relative strength 
of the nudging coefficient in this region. The assimilation has removed the eastward flow 
at 4°S, adjusting the flow in the T run (Figure 21a) towards the C run (Figure 21b). 



T C: U velocity componet 



SSH 156E (meters) 




Temp latitude 




-5 

^-10 

cd 

-15 

-20 



1.495-30.17 

Salinity 




-4-2 024 -4-2 024 -4-2 024 

latitude; 34.18 - 35.79 



Figure 21. Meridional section for 5°S to 5°N at 156°E of three model prognostic variables 
from the fields representing Sept. 5, 1992 T and C runs and the I field from December 28, 
1987. The top row (a-c) are the zonal velocities, the middle row( d-f) are the temperatures, 
and the bottom row (g-i) are the salinities. SSH is plot in (j), with the solid blank line from 
Sept. 5, 1992 of the C run, the dashed line from Sept. 5, 1992 of the T run, and the dotted 
line from the I field. The red line is the representative nudging term. 
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3. North Atlantic Ocean 



The statistical tests calculated for this region (Figures 9b-2, 1 Ob-2, 1 lb-2, and 12b- 
2) showed that the fields (T and S) are very similar. Plan views of the means of the 
variables T, S, U, V, and SSH also show little change (Figures 22 and 23, layers 1 and 
6). The December 92 map of standard deviations for run C (Figure 24 (layer 1) and 25 
(layer 6), row A) shows a much reduced variability of that over the initial fields (Figure 24, 
row C). With the assimilation, the variance is reduced (Figure 24, row B), unfortunately, 
the assimilation has over compensated for the reduction of SSH from the T run to the C 
run. It can only be speculated that in the adjustment process, the model needs a higher 
temporal sampling rate to include these smaller scale features. This is much more apparent 
in the surface layer (Figure 24) than at a depth of around 160m (Figure 25). 

4. Indian Ocean 

An interesting adjustment occurred in the assimilation run in the tropical region of 
the Indian Ocean. Plan views have been created of the various fields from the December 
1992 C and T runs and from December of 1988 of the C run, depicting dramatic 
differences (Figure 26). This figure shows the small area (20° x 20°) centered at 5°S, 90°E 
in the Indian ocean, with the top row. A, from the assimilation run, T; the middle row, B, 
from the control run C; and the last row, D, showing the fields from Dec. of 1988. This 
last row is provided to show the state of the simulated ocean, if starting from the same field 
which initialized the T run, but using winds from 1988, instead of 1992. This area 
corresponds to one of the boxes on Figure 6a and b, displaying the vertical nudging 
functionals. Because the wind forcing (Figure 26j and 26o) is not the same between 1988 
and 1992, a direct comparison can not be made as to how much of the change seen in 
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Figure 22. Plan views of monthly averaged u,v,T and S at 12.5 meters and SSH for the 
North Atlantic. The A row is from the Dec. 1992 C field, the B row is from the Dec. 1992 
T field, and row C is the initial or I field (Dec. 28, 1987). The range of the scale for each 
variable, u,v,T,S, SSH is listed below the image of row A. 
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Initial Dec. 87 




Figure 23. Same as 22, except for layer 6 at 160 meters. 
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Layer 6Depth: 160 





Figure 24. Same as 22, except for standard deviations. 
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Initial Dec. 87 




Figure 25. Same as 23, except for standard deviations. 
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Layer 6Depth: 160 






Figure 26- row A, is due to assimilation and how much is due to the difference in wind 
forcing. It is thought that a large portion of the change at the surface is due to the wind 
forcing because the nudging value for the top layer temperature in this region is relatively 
small. But changes are also seen at depth in this area (Figure 27, similar in layout as in 
Figure 26, but at a depth of about 600m), and the comparison of the fields of the control 
run for the same month (Figure 26, row B) are quite good in this region. In the 
temperature field, note the intrusion of a large amount of cold water centered at -5°S, 85 °E 
in Figure 26, row A, corresponding directly to a similar field in the control run (Figure 26, 
row B). Associated with this change in the temperature field are distinct changes in both 
the u and v velocity fields. At 600m (Figure 27), similar adjustments have been made in 
the T run (row A) to adjust the fields towards the control, run C (row B). 

D. CONCLUSIONS 

The twin experiment has shown that the assimilation process is adjusting fields both 
at the surface (SSH) and subsurface (temperature and salinity) towards the control 
simulation. Comparisons calculated from indirect estimates, such as those shown for 
potential vorticity, also lead to the conclusion that the model, if somewhat slowly, is 
adjusting itself towards the control. Four things have been brought to light in the twin 
experiment: 1) the T run has moved the model simulated fields towards what is observed in 
the control run and the next step is to run an experiment with real data; 2) most likely, there 
is a need for a stronger nudging term in the vertical in some regions; 3) higher temporal 
sampling may allow for better inclusion of mesoscale features; and 4) a longer period is 
needed to see the effect of the assimilation process when the nudging coefficient is as 
defined in this experiment. 
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Figure 26. Same as 22, except represent a 20’x20° box centered at 5’S, 90°E in the Indian 
Ocean. Additional column on right represents in row A, assimilation data; row B, wind 
field from Nov 92; and row C, wind field for I field, Nov. 88. 
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Layer 1 Depth: 12.5 










Figure 27. Same 26, except for layer 10 at 610 meters. 
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IV. ASSIMILATION OF ERS-1 AND T/P ALTIMETER 



DATA 

One goal in assimilating any observational data into a model is to correct the model 
towards what has actually occured. Then, once the model is correctly simulating the past 
and present, it can be analyzed as a 4-D dataset or be used to predict, with some accuracy, 
what the future state of a system may be. This has been the basic mode of operation for 
many years with weather forecasting. Although research in the area of assimilating 
observational data into ocean models has been active for twenty years or so, it is only 
recently that global data sets for the ocean with adequate temporal and spatial resolution for 
assimilation have been available. These data sets, for the most part, are created from 
measurements made from space, either of sea surface height (SSH) using altimeters or sea 
surface temperature (SST) using infrared instruments. For much smaller regions, 
instruments have been placed in the ocean for making long term measurements of the state 
of the ocean which are then used to adjust a model's fields towards reality. A good 
example of this type of experiment is the Tropical Atmosphere-Ocean (TAO) array set up in 
the tropical Pacific as part of the Tropical Ocean Global Atmosphere (TOGA) program 
{Hayes, et al. 1991). Extensive time series of current velocity, temperature and salinity 
have been collected for this data set and various researchers (e.g. Carton et al. 1996) have 
used the data to adjust tropical ocean models towards reality. 

On the global scale, the only complete data sets are those gathered from satellites. 
And, of the satellite data that can be used to adjust models, the best and most complete are 
those gathered by the altimeters, Topex/Poseidon (T/P) and ERS-1 and 2. SST 
measurements can not often be made because of cloud cover, thereby producing a 
somewhat patchy data set. There exists no global data set of temporally varying subsurface 
measurements on the scale of the POCM model, only climatological sets which are created 
by averaging observations and interpolating observations over many years. Thus, this 
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initial experiment of attempting to assimilate observational data into a high resolution global 
model is limited to the use of satellite altimeter data gathered from both ERS-1 and T/P 
satellites and subsurface information from climatological {Levitus 94) data sets. 

The rest of this thesis explains the assimilation process used with observational data 
and discusses the results of the experiment in light of three things. The first is whether the 
model's ocean circulation pattern is closer to the observations. The evaluation includes the 
use of previously described statistical tests to identify regions which have changed the 
most. In addition, dynamical balances are examined to determine if the assimilation 
process has changed the basic model dynamical properties. Finally, a comparison is made 
between the model fields from the control run (mn C), with no assimilation, and from the 
model run with assimilation (run A) to determine if the simulated fields have adjusted to 
observed fields in time and space. 

Second, the question: whether assimilating both ERS-1 and T/P data or T/P data, 
alone, are equally effective in pushing the model towards observations; is examined in 
Chapter V. The two satellites have different sampling patterns both in time and space. 
ERS-1 repeats its ground track every 35 days while T/P samples the same bit of ocean 
every 10 days. The implications of these two sampling patterns will be examined by 
looking at differences in the statistics of the model's variables of SSH, temperature, and 
salinity. 

Third and last, in Chapter VI, the adjustments to the model through the use of an 
assimilation process are examined in an attempt to understand the underlying deficiencies in 
the model simulation without any assimilation. For example, can the model reproduce the 
observed ocean with just a change in its initial conditions or are there some regions which 
will always be unpredictable? Is it important to simulate these turbulent features (such as 
Gulf Stream eddies) correctly or are statistically consistent fields over some temporal and 
spatial period accurate enough? This understanding will give a focus to why assimilation 
should be performed or not in attempting to answer the question: "what additional 
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information is gained about the global ocean's circulation that was not known before ?". 

The basic assimilation experiment is described next. Discussions of the quantitative 
and qualitative results of the assimilation of ERS-1 and T/P data into a high resolution 
primitive equation global ocean model complete this chapter. The chapter is then followed 
by two chapters which discuss items two and three above. 

A. SETUP FOR THE ASSIMILATION OF ERS-1 AND T/P 
DATA 

1. Description of the Assimilation Experiment 

The second experiment conducted for this thesis consisted of assimilating actual 
observed SSH anomalies as measured by two altimeters, the European ERS-1 satellite and 
the us/French T/P satellite. The experiment consisted of running the POCM, 
Semmer/Chervin model, with the inclusion of assimilation code, for the period from mid- 
1992 through the end of 1994. ERS-1 and T/P data were used jointly as observational data 
up to the end of 1993. After this time, ERS-1 was placed in an orbit with a temporal 
sampling period (168 day repeat), not conducive for use as an effective assimilation data set 
and only T/P is used for the observational field. A control run exists for the same period in 
which no assimilation has been performed, known as the POCM_4B model run and 
described in Stammer et al. (1996). The surface forcing used by both the POCM_4B run 
(C) and for the assimilation run (A) is the same: realistic 10m ECMWF wind stresses for 
the period 1992-1994 and monthly heat flux climatology described in Barnier et al. (1994). 
Both model mns, with and without assimilation, were performed at the National Center for 
Atmospheric Research (NCAR) on a Cray Y-MP, 8 processor computer. Although the 
assimilation period lasted from 1992 through 1994, the evaluation of the model output uses 
only the two year period spanning 1993 through 1994, thus eliminating some of the 
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adjustments the model is making initially to the model fields. 



2 . Handling of Altimeter Data 

The altimeter data for the assimilation experiment is from the French AVISO data 
set (Lc Traon et al. 1995) which was corrected for errors in the estimate of the ERS-1 
satellite orbit through the use of rnininiizing the distance between points (in time and space) 
sampled both by T/P and ERS-1 (i.e. cross-over analysis). These data cover the period 
from Oct. 6, 1992 through December 23, 1993 (T/P cycles 6-18). At that time, ERS-1 was 
in a 35 day repeat orbit and T/P has an established orbit of a 10 day repeat cycle. 
Atmospheric conditions have been corrected for in both data sets as described in Ch^ter I; 
and the tidal signal has been removed through the use of the empirical model, CSR 3.0 
(Eanes and Bettadpur, 1995). 

After the corrections where made to the along track altimeter data and a mean over 
the time period for the instrument was removed, the data were separated into files 
containing SSH anomalies for each day. These daily files, containing both ERS-1 and T/P 
data were then interpolated to the model's 1/4° mercator grid. If two satellite samplings 
occurred at the same point for the same day, the average of the two points was used. 

The earth's geoid is not known to the spatial resolution of the model, thus the 
altimeter SSH measurements are only the deviations of SSH from the 2 year temporal mean 
of SSH at a point (known as the SSH anomaly). But in the assimilation process, the total 
SSH signal is needed for the adjustment of the model to an observation. Thus, a SSH 
mean field needs to be added to the altimeter SSH anomaly values. Several different means 
can be used. One is a mean that is determined for in situ hydrographic data and another 
one would be an estimate of mean sea level firom altimetry itself (e.g. Rapp et al. 1996). A 
hydrographic mean is not easily understood because of the assumptions that have been 
used for determining a level of no motion. The use of an estimate made from altimetry is 
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good only to scales on the order of 2000 km, much larger than the model's grid size. Rapp 
et al. showed that the altimetric mean at these long wave lengths compared well to the 
POCM_4B model, therefore the SSH mean from a previous run of the model (POCM_4A, 
an average of three years, 1987 through 1989) will be used as an estimate of the mean 
SSH. The use of a mean which is not the mean of the control run was done to examine 
how a somewhat different mean SSH field is handled by the assimilation process. A test 
model run which used the mean SSH field from the control mn was also conducted for 
1993. The analysis is contained in the appendix. The two mean fields have a global 
difference in their mean of 1.9 cm +/- 7.8 cm. The mean SSH (<SSH>) field added to the 
altimetric values is shown in Figure 28a with the mean SSH from the C run for years 1992 
through 1994 (corresponding to the assimilation period) shown in Figure 28b. 

The assimilation process updates the altimeter observation field each day. If a new 
observation appears at a grid point, it is given an associated temporal weight of one. For 
all the observations which are not updated from the newly read daily observations file, the 
weighting factor associated with an observation is increased by one day, up to a maximum 
of 30 days. 

3. Associated Vertical Relationships and Nudging Coefficients 

Next, SSH observations need to be associated with the vertical structure of 
temperature and salinity. The temperature and salinity functions which relate SSH to 
subsurface values (Equations 2.5 and 2.6) are created as described below for this second 
experiment. Monthly maps of joint ERS-1 and T/P SSH fields (including the SSH model 
mean) are created from the joint AVISO data set A 5x5 point smoother is then applied to 
each monthly SSH field to remove very high wave number features which do not appear in 
the temperature and salinity data set used for forming the functionals. The subsurface 
fields of temperature and salinity used are from the monthly climatological fields of 
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Figure 28. Mean SSH fields for a) 2-year (93/94) C run and b) 3-year (87-89) P0CM_4A 
used for SSH mean to add to altimetric SSH anomaly observations. 
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temperature and salinity as produced by Levitus and Boyer (1994) and Levitus et al. 
(1994). The Levitus fields are improved from the 1982 data {Levitus, 1982) and should be 
somewhat dynamically consistent with the observational altimetric data set, since both 
result purely from measurements of the real ocean. The Levitus 94 data set is not the data 
set which initialized the ancestor of the current POCM_4B model (a one-half degree model, 
see Semtner and Chervin, 1992), nor is it used to keep the model close to the observed 
values in the boundary regions. Thus, it provides a set of observed temperatures and 
salinities (albeit averaged over different time and space periods ) which have not been 
incorporated into the model previously. The other option available would have been to use 
fields of temperature and salinity from a previous model run or the control period (1992- 
1994 of POCM_4B). It is thought that if a set of temperature and salinity data from a 
previous model run were used to form these functionals, the assimilation process could not 
be fully tested because the vertical structure would be too similar to the model run without 
assimilation. 

Covariances computed between the monthly SSH fields and the monthly 
temperature and salinity fields as defined by Equations 2.5 and 2.6, respectively, are then 
normalized by the variance of the SSH at each grid point. Figures 29a and b are plots of 
the functionals for the second experiment similar to the Figure 6 for the first experiment, 
the twin experiment. As is seen, the functions are a somewhat smoother than the functions 
created for the twin experiment which used 3-day fields over a year to create the 
covariances. 

The last field to be constructed for use in the nudging process is the field that is 
related to the e term in Equation 2.2. This term is the normalized variance of SSH 
computed from observational satellite values gridded into two degree boxes and with a 5x5 
grid smoother applied. The variances are normalized by 40 cm^ (chosen empirically), with 
any value over 1 , set to 1 . 40 cm^ is close to the maximum variance, but allows for the 
elimination of the extremely high variances that occur in the Irish Sea and English Channel. 
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Figure 30 shows the spatial structure of the weighting field. The model is pushed very 
aggressively towards observed values in the region of the Kuroshio, Gulf Stream, and in 
parts of the Antarctic Circumpolar Current. This map could be made to reflect a much 
different weighting. For example, it could be all one value or created from simulated fields 
instead of observed fields. Keeping in mind that the assimilation with observed SSH fields 
over the globe at high resolution is the first to be attempted, an initial guess to define this 
relaxation weight field is an that is warranted at this time. Part of the evaluation process is 
to determine if this field should be modified. 

Several model runs of the assimUation process were performed initiaUy to produce a 
realistic simulation for evaluation. The underlying nudging coefficient (see Equation 2.2) 
basic relaxation coefficient, RO, was given several values. (In reaHty, e = RO times the 
normalized variance of SSH, where RO is some time reciprocal.) Determined, empiricaUy, 
the RO nudging coefficient of the temperature and salinity fields is equal to 30 days. This is 
directly related to having used monthly fields for the creation of the vertical fields. For the 
surface fields, RO is set at a quarter of a day. 

Three lines are shown on Figure 31a representing the total model energy A run 
(black line), the C mn (blue line) and an earlier run with RO equal to half a day (red line), 
while Figure 31b and c show the total rate of change in the temperature and salinity fields. 
Only when RO is increased to 0.25 days (at time = 1993.75 on the red line), does the 
energy in the model increase. Because of the changes shown between the blue and red line 
near the end of 1993, the final assimilation run, run A, used RO equal to one quarter of a 
day for height. It should be remembered that this value is only for very high variable 
regions. Regions which show much smaller variances in Figure 30 are not nudged nearly 
as strong. 

Two final checks had to be included in the assimilation process to limit the 
maximum allowable change in temperature and salinity. The problems exhibited 
themselves when the model became unstable at depth. The change in the observed and 
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modeled temperature at a grid point is, therefore, not allowed to change by more than 5°C. 




Figure 30. Spatial structure of weighting field for nudging coefficient. Scale is from 0 to 1 . 



71 



ENERGY 




Figure 31. Time series of diagnostic quantities during model run for a) total energy, b) 
change in temperature, and c) change in salinity. Black line is the A run; blue defines the 
control run, and red line are the values from an early run with differing nudging 
parameters. 
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B. STATISTICAL EVALUATION 



1. Comparison of Means and Variances 

The statistical tests used for evaluation in the twin experiment are used to examine 
the magnitude of the changes in the second experiment, involving the assimilation of real 
altimeter data into an ocean model. Simple comparisons of the differences in the two year 
(1993/1994) means between the A and C runs of SSH, temperature, and salinity quickly 
elucidate where the model has changed the greatest. At this point, nothing is being 
indicated about whether the changes produce a more realistic model. The discussion 
comparing the model runs to observations is in section D of this chapter. The differences 
greater than 10 cm between the means of SSH for the two runs are shown in Figure 32a. 
Changes are clearly seen in regions where they are expected to be seen. In the Kuroshio 
Extension of the North Pacific, the assimilation of the altimeter data has increased the mean 
sea level, while in other areas, such as the south subtropical Indian Ocean, the mean has 
decreased. Throughout the rest of the domain, the significant changes alternate between 
regions where the C run is higher than the A run and vice versa. Figure 32b shows 
differences between the C and A runs for the standard deviation of SSH with only 
differences greater than 5 cm outlined. This figure clearly shows that with the assimilation 
of altimeter data, the SSH field has become more variable in the Southern Ocean, the 
Kuroshio Extension, and in the Gulf Stream region. The one region which is showing a 
decrease, is the path which the Aghulas eddies take in the C run of the model centered near 
35“S 10°E. The model C run shows a more northerly path of high variability over which 
eddies can be tracked. In the A run, the path of the eddies (as defined by the altimeter data 
being assimilated) as not so distinct. The variability pattern (here represented by the red 
contours centered at 30°S 10°E) is much as is seen in the variability maps calculated from 
altimetry data alone. The differences in the eddy kinetic energy (EKE), defined as; 
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[(u-u)2+(v-v)2], with the overbar denoting the mean value, for the top layer between the A 
run and C run is shown in Figure 32c. The EKE field, computed from the velocity 
components, is a measure of the how the assimilation of SSH is affecting the variables of 
the model which are only indirectly, through the model's physical equations, related to the 
variables which are being modified through the assimilation process. Similar to the 
changes seen in the SSH fields, the EKE has been modified in the regions which are 
associated with high variability, the Southern Ocean, the Gulf Stream, the Kuroshio 
Extension, and the tropics. This can be traced to equation 2.2, in which part of the 
relaxation coefficient, e, is a function of the variability. 

Further down in the water column, at 160 m or level 6, changes in the temperature 
field can also be seen in both the mean field and the standard deviation field (Figure 33a,b). 
And again, the changes are the greatest where it is expected, in the regions of high SSH 
variability. 

2. Further Statistical Tests 

As calculated for the first experiment (twin), statistical tests are used to quantify the 
changes seen between the C and A run of the model. Figure 34 shows the student-t and F 
test values calculated between various pairs of monthly fields, globally and for the same 
regions as shown in Chapter 11 IB-2. The thin lines give the values computed between the 
January 1993 (year 1) temperature and salinity fields for the C and A runs. The bold lines 
are for January 1994 fields (mid- way point) and the dashed lines are at the end of the 
experiment December 1994. These comparisons have been calculated using monthly mean 
fields, rather than the snapshot fields shown for the twin experiment. The horizontal axis 
for Figure 34a and b is defined to have the SSH t and F test values at x = 1, with the 
temperature layers 1-19 shown at x=2-20. It be should expected that the greatest difference 
should be in December 1994. The t-test for temperature and SSH, reflecting how similar 
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Figure 32. Differences in the SSH fields between the A and C run for a) the 2 year mean 
(93/94) and b) the 2 year standard deviation. Red defines areas where the A run is greater 
and B defines areas where the C run is greater. 
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using SSH. Red defines areas where the A run is greater and B defines areas where the C 
run is greater. 
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Figure 33. Differences in the temperature field at 160 meters (layer 6) for a) the 2 year 
mean field and b) the 2 year standard deviation field. Red defines areas where the A run is 
greater and B defines areas where the C run is greater. 
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or different the fields are in the mean, shows the global mean (Figure 34a- 1) is statistically 
different only near the surface. Some of the regions, though, have more significant 
changes in the mean, for example in the Indian Ocean (Figure 34a-6) and the tropics 




Figure 34a. Plots showing the Student-t test values calculated between different model 
fields for the globe and for different basins for temperature. A pair of panels comprise each 
figure. Y-axis is the t-test value, 1.64 is the significance level for a null hypothesis. X-axis 
is read as follows: x = 1, t-test value for SSH, x=2-21, t-test value for temperature, levels 
1-20. There are six areas: 1) global values 2) North Atlantic, 3) South Atlantic, 4) North 
Pacific, 5) South Pacific, g) Indian, 7) Tropics, 8) Southern Ocean. The t-test values are 
computed between a pair of model fields, one from the A run and one from the C run. The 
bold line represents the month of August 1993 (midway in the assimilation run); the dash- 
dot line represents the month of December 1994 (end of the assimilation run); and thin line 
represents January 1993 (early in the assimilation run). 
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(Figure 34a-7). This corresponds directly to the changes seen in the plan view of the 
differences shown in Figures 32a (surface) and 33a (160 m, level 6). Even in the North 
Pacific, the t-test value for level 6 (x=7) in Figure 34a-4, shows an increase in the t-test 
value for the final month (dashed line) of the experiment. The F-test calculations. Figure 
34b, shows that the variability is significantly different between the two model runs, with a 
larger change occurring during the first year than in the second (Figure 34b- 1). The North 
Atlantic is the only region that doesn't show the December 1994 change (dashed line) 




Figure 34b. Same as 34a, but are the F-test values for temperature. 
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greater than the December 1993 changes (bold line) for top 1500 m (level 12). below this 
level, it is difficult to understand the significance of the t and F test values because the 
changes, when they do occur, are very small. 

Similar computations have also been done for the salinity fields for the same 
monthly pairs (Figure 35a,b). The changes do not seem to be as significant for salinity as 
they were for the temperature fields, although over some regions, such as the Southern 
Ocean (Figure 35a-8) and the Indian Ocean (Figure 35a-6), there does seem to be 
significant change in the mean of the salinity field. The F-test calculations (Figure 35b) 
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In summary, this section has shown, through the use of statistics, that the model 
has changed significantly both at the surface and at depth, with the inclusion of 
observational data via an assimilation process. Further, regional analysis has shown that 
these modifications to the model's prognostic variables are not the same for all regions and 
the change is most directly related to the variance in the observed SSH field. 
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C. DYNAMICAL CONSIDERATIONS 



For an assimilation process to be successful, it should not change the basic 
dynamics as defined by the model's physical equations. Two things are examined to 
ascertain whether the inclusion of the assimilated SSH fields and associated vertical 
structures, has changed the dynamical properties of the model. The balances between the 
various forces are looked at first, followed by a section which determines if potential 
vorticity along several Montgomery streamlines {Le. streamlines defined on an isopycnal 
surface) is consistent between the two runs of the model. 

1. Balance Between Forces 



The u-component of the momentum equation is written as 



3u ^ 3u 3u 3u ^ ^ 1 3P 3^u . * /3^u 3^u>^ 



4.1, 



where k and A are specified fiiction coefficients. The two year mean fields of the A and C 

runs are used to determine the order of magnitude of the above terms to determine if the 

assimilation process has modified any of the physical balances defined by the model's 
numerical equations. The advection term, u^ *he Coriolis term, fv, the pressure 

gradient term written as a function of SSH (via the hydrostatic approximation), gg^ , and 
the horizontal friction term, are computed separately and then compared. The 

calculations were performed for the whole globe and can be represented using two 
meridional bands averaged over 10*, one in the Pacific at 157°E and one in the Atlantic at 
32.5°W. Figure 36b shows the magnitude of the each of the four terms centered at 
157.5°E, with the dashed line representing the A run and the solid line, the C run. From the 
plots, it is difficult to see any diversions of the dashed line from the solid line. Thus, the 
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model output from the two runs give approximately the same values for all the terms (as 
expected from scale analysis) except for the pressure gradient and Coriolis terms near 




Figure 36. Zonal average over a 10°xl0° box of the pressure gradient, Coriolis, advection, 
and friction terms of the momentum equation, centered at a) 327.5 °E in the Atlantic and b) 
157.5°E in the Pacific. The solid line is the C run; the dashed line is the A run. 
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45°N and 45°S and in the Southern Ocean for the friction term. Most important, though, 
the order of the terms do not differ between the two mns of the model. The meridian in the 
Atlantic at 327.5°E shows a similar comparison (Figure 36a). 

2. Potential Vorticity 

Potential vorticity (PV) can be used as a diagnostic quantity to evaluate the changes 
to the dynamics of the model when SSH and a related temperature and salinity structure are 
assimilated. Gille (1997), determined, using model output from the Semtner/Chervin 
POCM model, that PV is not conserved along Montgomery streamlines in the Southern 
Ocean. Even with this knowledge of PV behavior along a streamline, potential vorticity 
can be used to determine if the two runs of the model are dynamically similar by looking 
for differences in the order of variation in PV between the mns. 

An isopycnal surface is a surface which is more dynamically suited for looking at 
flow properties because the flow is maximized along a surface and minimized across 
surfaces. The Montgomery streamfunction {Montgomery, 1937), the geostrophic flow 
along an isopycnal and defined as a pressure anomaly term by Zhang and Hogg (1992), 
can be used to represent the geostrophic relation for the horizontal geostrophic velocity. 
The constraint of a non-divergent flow does not hold as it does in the traditional definition 

A 1 

of a streamline. The geostrophic relationship is defined d&k \fn =—Vzp. This equation 

P 

can be transformed into an equation describing geostrophic flow on an isopycnal by 
substituting in the hydrostatic equation (^ = -pg) and using the equation, V 2 A = V 0 A 

VqA, to transform the geostrophic equation to isopycnal coordinates, resulting in: 

^ x/ii =- VoP - Vaf-dp 4.2. 

P P 

Zhang and Hogg analyzed the equation and determined that the non-gradient component 
can introduce significant error in calculating geostrophic velocities along an isopycnal, 
when the total pressure is used. If the mean pressure is removed from the pressure term. 
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the resulting equation substantially reduces the error and gives a more reasonable estimate 
of the velocity field on an isopycnal. Therefore, in this paper, the Montgomery 
streamfunction is defined using the pressure anomaly instead of the total pressure with the 
equation written as: 

M = 5(z) (p(z) - p) - 5(z) dp 4.3 , 

where z represents the depth of the isopycnal surface and 6(z) is the volume anomaly 
defined as 



6(z) = 



1 

P(z) 



1 

Po(z) 



4.4. 



p(z) and Po(z) are the in situ density and the reference density with salinity equal to 34.6 
and temperature equal to 2.33°C. For these calculations, the mean pressure removed is 
equivalent to around 300 meters. The free surface is included by modifying 4.3 to 



M = SpsTig + 5(z) (p(z) - p) 



- F 5(z) 



dp 



4.5, 



with ps equal to the surface density and ti is the surface height. 

do 



PV is calculated as Q = f — from the model's prognostic temperature and salinity 

dz 



values and is calculated at the points defined by the Montgomery streamline. Relative 
vorticity has been eliminated from the computation because GiUe (1997) has shown that it 
contributes little to total PV, most of which is determined by the vertical structure. 

To examine if the assimilation mn is dynamically consistent with the control run, Q, 
plotted along Montgomery streamlines, is examined for several areas. The flow is 
examined only to determine if the assimilation values of potential vorticity are consistent 
with the control run. At a depth of around 1300 meters. Figure 37a shows the 0.09 
Montgomery streamline in the Southern Ocean as calculated from equation 4.5 from the 
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Figure 37. Various quantities along a Montgomery streamline on a density surface of Oi = 
32.3 in the Southern Ocean, a) shows the streamline path for the A run (solid line) and the 
C run (dashed line); b) is the potential vorticity along the path; c) the thickness of the 
density layer, and d) the depth of the layer. 
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mean of 1994 along the 32.1 Oi surface. The heavy line represents the quantity as 
calculated from the C run and the thin lines are those computed from the A run. Panels b-d 
show the PV, the thickness and the depth of the layer, respectively. Changes have 
occurred in the model run at this depth directly corresponding, in most cases, to the 
thickness of the layer, especially between 140° and 150°E, near Tasmania. The minimas 
and maximas of both runs are of the same order. Around 10°E, the streamline has changed 
dramatically, but there is little change in the PV along the line. This is the path of the 
eddies shed by the Agulhas Current into the South Atlantic. This evidence suggests the 
two runs are producing different, but dynamically similar flow at depth. In other words, 
the assimilation is not modifying the underlying dynamics of the model. 

As we look further up in the water column in the Southern Ocean, a specific 
streamline shows more differences in location between each of the model runs. For 
example. Figure 38a-d shows the changes occurring along the -0.05 contour on the oi 
surface in the Southern Ocean. The magnitude of the PV along the line remains essentially 
the same between the two runs of the model. At a depth between 600 and 800 meters, the 
-0.13 streamline has also been modified between the two model runs, with much larger 
changes as seen in Figure 39a-d. Above this depth, it becomes increasing more difficult to 
computationally follow a streamline. This is also true for other regions, such as the North 
Atlantic and the Kuroshio area of the North Pacific. 

In summary, these d 5 mamical comparisons of the balances in the terms of the 
momentum equation and of the consistency in the flow along streamlines in the Southern 
Ocean show that the assimilation process is not modifying the basic physics of the model's 
numerical equations. 
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Figure 38. Same as 37, but for oi = 32.0 in the Southern Ocean. 
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Figure 39. Same as 37, but for Oi = 31.8 in the Southern Ocean. 
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D. COMPARISONS TO OBSERVATIONS 



As stated above, one of the goals of assimilating observations into a model is to 
adjust it towards actual observations. This section looks at several things to determine if 
this goal is being met. The sea surface is examined first by comparing several quantities 
with observations. This includes examining the differences in the variability and spectra of 
the surface fields as compared to altimetric data. Also included, is the comparison of time 
series of model SSH at specific locations to observational tide gauge measurements. 
Subsurface, the changes in the heat content of the model are examined, along with looking 
for changes in the thermohaline circulation pattern when observational data is assimilated. 
This is in addition to comparing the model's diagnostic fields (velocity, temperature, and 
salinity) with observations sampled at the same time in the ocean, using observations 
mainly retrieved from the World Ocean Circulation Experiment (WOCE) database. 
Because the domain of the model is global, this thesis will concentrate on examining the 
changes in the regions where the assimilation is expected (using the spatial structure of the 
weighting function. Figure 30 as a guide) to have modified the dynamics the most. These 
areas are 1) the North Atlantic in the region of the Gulf Stream, 2) the Kuroshio Region in 
the North Pacific, 3) the tropical ocean, and 4) the Southern Ocean. Results from the 
evaluation are shown from other areas when they are pertinent to the discussion. 

1. Sea Surface Height 

Figure 32a shows the gross differences in SSH between the two year means of the 
two model mns, C and A. The SSH field can be examined in closer detail by focusing on 
a region. The four regions mentioned above are examined to determine if the estimate of 
SSH has improved in the A run (the assimilation run) over that of the C run (the control). 
A time series of the RMS difference between the SSH in mn A and the observational SSH 



90 



(Figure 40), shows that over the first 150 days of the model run, the RMS difference in the 
observed field drops from over 16 cm to about 13 cm, where it remains for the remainder 
of the two year run. (If the mean SSH from the C run is used as the mean SSH to add to 
the SSH anomaly field, the RMS difference drops only slightly, to 1 2 cm.) A large part of 
the difference remaining after day 150 is related to the regions of the model which have 
high variability that is not seen in the observations, but appears in the model fields. 




Figure 40. RMS difference between model SSH and the observations from day 0 of A run. 



Because the assimilation process is only adjusting those areas where there is high 
variability in the observations, these other high variable areas largely remain unchanged. 
This is, perhaps, where a correction should be made to the assimilation process. If the 
nudging coefficient R, as used in Equation 2. 1 , were a constant equal to the maximum of 
the current R, instead of a value dependent on the variance of the SSH at the location, the 
assimilation process should adjust the model towards the model mean faster in these areas 
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that have changed relatively little with the present assimilation scheme, this is because the 
nudging weight is related to the variance of the SSH, not the mean field. The current R 
value allows for changes with higher variance to be nudging more strongly or in other 
words, with a shorter relaxation time scale. At this point in time, it is difficult to say which 
is preferable, partly because the mean used for the observational mean is a model mean. 
The next paragraph looks at the adjustments to the mean and standard deviations for the 
regions mentioned above. 

Figures 41 through 43 show the standard deviation of SSH for the control field (a), 
for the observation field (b), and for the assimilation field (c) for 3 of the 4 regions. (The 
tropical region shows little variance change and is not shown.) The Kuroshio region. 
Figure 41c, shows clear enhancement of the SSH anomaly field in the latitude band 
between 30°N and 40'’N matching the variability seen in the observations (Figure 41b) with 
the highest variance. There is still a diagonal band of high variance, not seen in the 
observations, from about 145°E, 40°N to 160’E, 45°N. In the Gulf Stream region, the 
band of high variability is flattening out somewhat in the A run (Figure 42c), but not as 
dramatically as it should, as seen in the observed field (Figure 42b). The Southern Ocean 
(Figure 43a-c) shows similar variance stmcture in those regions where the model has been 
nudged the strongest (40’S, 0° and 70’E, 60“S, 180’E). The assimilation fails to duplicate 
the variance structure in the Brazil-Malvinas Confluence (45°S, 310°E). 

The means of these regions are shown in Figures 44 through 46. The mean of the 
Southern Ocean has changed little and therefore, is not shown. The Gulf Stream region's 
mean field from run A shows the permanent eddy-like feature located in Figure 44a at 35 °N 
289°E as either much reduced or removed completely. SSH plotted along the 290°E 
meridian shows this change as a change in the location of the Gulf Stream wall (Figure 
44d). In the Kuroshio region, the mean contour of -50 cm centered at 50°N, 165°E is 
enlarged in area in the assimilation (Figure 45c) over the contour in the control run (Figure 
45a). This suggests that the assimilation mn is closer to the observed mean field (Figure 
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Figure 41. Contours of the 2 year standard deviation of the SSH field for the Kuroshio 
Region; a) run C, b) observations, and c) run A. . Contours are every' 5 cm. 
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Figure 42. Same as 41 , except for the Gulf stream region. 




Figure 43. Same as 41 , except for the Southern Ocean 
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Figure 44. Mean field from the Gulf Stream region. Contours are 50cm. a) the C run, b) 
the observations, and c) the A run d) shows the mean SSH along 290°E for the C run 
(dashed line) and the A run (bold line). 
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Figure 45. Same as Figure 44a-c, but for the Kuroshio region. 
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Figure 46. Same as Figure 44a-c, but for the global tropical regions. Contours are 40cm. 
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45b) than is the SSH from run C. The contours of 0 and 50 (between 35° and 40°N) 
appear to be closer in the A run to the observed contours than is seen in the C run. In the 
tropical region (Figure 46), the 50 cm contour has moved westward in the Pacific ocean in 
the A run from where it is located in the C run. The A run contour (Figure 46c) is closer 
to the observed field (Figure 46b) than is the C run (Figure 46a). These plots show that 
the mean of the model has changed along with the variance of it, adjusting it towards the 
observed mean. 

Sea level heights, as measured by tide-gauges, have not been assimilated into the 
model, and thus, can provide an independent measure of the improvement in the model 
with the assimilation of altimetric SSH anomalies. Before showing the results of the 
comparisons of the assimilation run to the tide-gauge measurements, the performance of the 
control run against these independent measurements (when no observations have been used 
to adjust the model's prognostic variables) is examined to show that even without 
assimilating observations, the sea level of the model is simulated adequately. 

Tokmakian (1996) gives an in-depth comparison of POCM_4B SSH values to sea 
level measurements of 94 stations over a period of 9 years. Time series of the model SSH, 
averaged over 60 km and sampled every three days, were created at the locations of the 
tide gauges. Similarly, the sea level data from the tide-gauges were sub-sampled from the 
daily values (with a seven day low pass filter applied) to a 3-day time series to correspond 
to the same sampling of the model. An inverse barometer effect was removed also from the 
observational data. Finally, the two year mean of each data set was removed before the 
correlations were computed. The overall mean correlation is 0.48 +!- 0.28 for 1993- 
1994. At some of the stations, the analysis of the C run shows that the model is 
reproducing the observed signal quite well. For example, the plots of the time series at 
Naos Island (Gulf of Panama, 8.9°N 79.5°E) and at Christmas Island (2°N 157.5°E) 
(Figure 1) shows the El Nino events of November 1991 and 1994 at Christmas Island 
which are then seen about a month later at Naos Island. This fits the estimated wave speed 
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of (2.8 m/s) given be Delcroix et al. (1994). 

To evaluate the improvement to the model with the assimilation of the altimetric 
data, the same comparison is done with run A. The resulting correlations (with an annual 
cycle removed from both observations and modeled SSH) to the tide-gauges have been 
divided into three sets 1) where there is little change in between the C and A runs of the 
model (to within 0.1 ), 2) where the A run performs better, and 3) where the C run 
performs better. There are 59 stations where there is little change between model runs. Of 
these stations, 50 have a correlation better than 0.3 and 41 are over 0.5. In the second 
group (Table 3), locations at which correlations were less than 0.15 were removed, 
leaving 19 locations with an improved correlation in the A run. The third group is shown 
in Table 4. Here 17 locations have worsened with the inclusion of the assimilation 



Tide Gauge 


R; A Run 


R: C run 


Latitude 


Longitude 


POHNPEI 


0.62 


0.44 


7.0'N 


158.2°E 


MAJURO 


0.70 


0.54 


7. UN 


177.4 


MALAKAL 


0.19 


0.03 


7.3'N 


134.5 


FRENCH FrSh 


0.36 


0.05 


28.9'N 


166.3'W 


RKITEA 


0.51 


0.34 


23. rs 


134.9°W 


SUVA 


0.71 


0.54 


i8.rs 


178.4'E 


ADAK ISLAND 


0.29 


-0.02 


51.9°N 


176.6°W 


CHICHUIMA 


0.46 


0.32 


27. UN 


142.2‘E 


MIDWAY ISLAND 


0.60 


0.11 


28.2'N 


177.4"W 


JOHNSTON 


0.33 


-0.01 


16.8“N 


169.5“W 


HILO 


0.66 


0.12 


19.7'N 


155. rw 


CHATHAM 


0.52 


0.18 


44.0^8 


176.6^W 


ST.PAUL 


0.73 


0.58 


38.7“S 


77.5‘E 


KERGUELEN 


0.56 


0.31 


49A°S 


70.2‘E 


KEY WEST 


0.31 


0.10 


24.6^N 


81.8W 


SETTLEMENT PT 


0.31 


0.14 


26.7“N 


79.0"W 


PORT STANLEY 


0.48 


0.27 


51.7°S 


57.8'W 


ABURATSU 


0.24 


0.03 


31.6'N 


131.4^E 


NAHA 


0.70 


-0.15 


26.2'N 


121, TE 



Table 3. Tide-gauge stations where the correlation, R, is higher for A run than C run. The 
correlations listed cover the 2 year period 1993-1994. An annual harmonic has been 
removed. 
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Tide Gauge 


R: A Run 


R: C run 


Latitude 


Longitude 


JUANBERNANDEZ 


0.45 


0.61 


33.6°S 


78.8®W 


PENHRYN 


0.30 


0.48 


9.0*’S 


158.rw 


HONOLULU 


0.17 


0.53 


21.3°N 


157.9‘W 


VALPARAISO 


0.45 


0.61 


33.0°S 


71.6°W 


PORT LOUIS 


0.06 


0.38 


20.2°S 


57.5°E 


RODRIGUES 


0.14 


0.58 


19.7°S 


63.4°E 


CALDERA 


0.27 


0.47 


27. rs 


70.8"W 


DARWIN 


0.28 


0.45 


12.5'S 


130.8°E 


ESPERANCE 


0.27 


0.68 


33.9°S 


121.9°E 


DUCK PIER 


0.40 


0.70 


36.2°N 


75.7®W 


OFUNATO 


0.07 


0.42 


39.1°N 


141.7°W 


CRESCENT CITY 


0.57 


0.67 


41.7°N 


1 124.2^W 


NEAH BAY 


0.56 


0.71 


48.4°N 


124.6°W 


SAN DIEGO 


0.53 


0.64 


32.7°N 


117.2“W 


YAKUTAT BAY 


0.55 


0.71 


59.6°N 


139.7'W 


KETCHIKAN 


0.44 


0.62 


55.3°N 


131.6'W 


DIEGO RAMIREZ 


0.25 


0.44 


56.5°S 


68.7°W 



Table 4. Tide-gauge stations where R is lower for A run than C run. The correlations listed 
cover the 2 year period 1993-1994. An annual harmonic has been removed. 

process. When the annual harmonic is included in the computation of the correlations 
between time series, eight of the stations in list 3 move to either list 1 or 2 (Darwin, 
Esperance, Duck Pier, Ofunato, Crescent City, Neah Bay, Yakutat Bay, and Ketchikan), 
indicating that the high frequency variability along coasts is being degraded by the 
assimilation process. None of the stations in list 2 move to list 3 if the annual harmonic is 
included, meaning the annual cycle is not degrading the correlation at the locations 
improved in run A. The annual harmonic does degrade the correlation values for the A run 
at three stations, Easter Island, Kushimoto, and Saipan which showed had shown no 
change in the correlation values for higher than annual frequency variation. 

By looking at individual locations, a qualitative evaluation can be made in how the 
models simulate the sea level. At Naha (26.2'N, 127.7°E), in the Kuroshio region, the 
correlation between the time series of tide-gauge data (with an annual harmonic removed) 
and the model increases from -0.15 to 0.70 (Figure 47a). At Naha, it is mostly the long 
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temporal scale (seasonal or greater) of sea level variation that has increased the correlation 
value, clearly seen qualitatively in the plot with run A (heavy line) appearing to follow the 
tide-gauge data (thin line) better than run C (the dotted line). In the tropical region at 
Pohnpei (7°N, 158.2°E, Figure 47b), the correlation between the C run and the tide-gauge 
measurements is 0.44 and with the assimilation of the altimetric data (run A), the 
correlation increases to 0.62. The increase at this location is mostly an improvement in the 
higher frequency (less than a season) events which occur in 1994, especially around the 
October/November time frame. In the middle of the North Pacific at Midway Island 
(Figure 47c), the assimilation improves the model from 0.11 to 0.60, with most of the 
improvement occurring at high frequency. Note here, the magnitude of the model SSH 
anomaly still underestimates the observed sea level anomaly. 
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R Assim 0.7 R Control:-0.15 

Ai 




Location: 26.22N 127.67E station: 89 





YEAR 



Figure 47. Times series of SSH for 3 stations: a) Naha in the North Pacific, b) Pohnpei in 
the tropical Pacific, and c) Midway Island in the North Pacific. The bold hne represents the 
A run; the thin line is the data measured with a tide gauge; and c) the dotted line is the SSH 
from the C run. An annual harmonic has been removed from all the time series. 
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2. Spectra of SSH, Surface Velocity Fields 



The differences between the SSH and velocity fields can be further quantified by 
examining the spectra of these fields. Average 10°xl0“ spatial SSH meridional spectra are 
computed as in Stammer and Boning (1992), and then averaged over a year from 10 day 
maps of SSH. A map of every other average spectra is shown in Figure 48 for the whole 
globe. From the global viewpoint, when averaged over these 20°x 20° areas, the spectra 
between the two model runs has changed little. If we look at several locations in particular, 
at the region centered around 45°N, 75°W in the Gulf Stream area and the Kuroshio region 
(45°N,185°E), the observational spectrum (dotted line) (Figure 49a,b) is different from 
either of the two model mns, but the stiU within the error of the observed spectra. At 35 °S, 
5°E (Figure 49d), the assimilation has reduced the variance of the model at low wave 
numbers, more representative of the altimetric observations. In the western tropical Pacific 
(Figure 49c), it appears that the assimilation has not changed the spectrum. 

The zonal velocity temporal spectra are calculated from 3 day fields and along two 
meridians, 30° W in the Atlantic and 200°E in the Pacific at approximately every 10 degrees 
of latitude. No averaging over latitude bands has been done for these spectra. Figures 50 
and 51 show the spectra as it is computed from the C run (thin line) and from the A run 
(bold line). The figures confirm what has been seen in the change of EKE of the model 
(Figure 32c). The A run shows a significant increase in the velocity variance over the 
whole spectra in the Southern Ocean (50°S and 30°S) in the Pacific. In the tropics at 10°S, 
the A run spectra is slightly higher, but not significantly. While in the North Pacific, there 
is little change in the spectra along this meridian. This is also consistent with Figure 30, 
showing the locations where the model is nudged the strongest. There is weak nudging 
along this line in the 200°E meridian in the Pacific. 
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Figure 48. Global spectra of SSH of the A run (black line) and C (red line) averaged over 
20°x 20° areas. The left side labels indicate the latitudes and the labels across the bottom 
indicate longitude. The y axis of each small plot ranges from 5e-4 to 5e-2 and the x axes 
range from le-1 to le6. 
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Figure 49. SSH spectra for the a) Kuroshio region, centered at 45 °N 185°E; b) for the Gulf 
Stream region, centered at 45°N 285°E; c) for the western tropical Pacific, 15°S 165°E; and 
d) for the Aghulas eddy path in the South Atlantic, centered at 35°S 5°E. The dashed line 
are the C run values, the dotted line is from the A run, and the solid line is the observed 
spectra. 
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Figure 50. Spectra of the zonal velocity along 200 ""E in the Pacific, at every lO"" of latitude. 
The bold line are the A run values and the thin line are from the C run. 
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Figure 51. Same as 50, except for the SSO^’E meridian in the Atlantic. 

3. Comparisons of Subsurface Prognostic Variables to Observations 

The four regions described above are examined again in this section, but now 
looking at the changes that have occurred in the vertical. The subsurface fields are 
compared to observational hydrographic data to determine if the vertical relationship 
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between SSH and temperature or salinity improves the simulation of the subsurface fields. 
Each area is examined in two ways, first by qualitatively and quantitatively describing the 
changes in the vertical structure and second by comparing the model sections to an in situ 
hydrographic section. 

a. Kuroshio Region 

Figures 52 and 53 show plan views of five prognostic model variables, 
U,V,T,S, and SSH (t)), along with the observation field for the month of December 1993 
(approximately half way through the assimilation run). SSH is included for completeness 
to show how well the observation SSH has been included into to run A's SSH field. 
Figure 52 shows the fields at a depth of 37.5m (level 2 of the model) while Figure 53 
shows the fields at a depth of 160m (level 6) and the surface height field (ti). The bottom 
panels are for the C run and the top panels are for the A run. Below the first set of panels 
is the color range for the values for a particular variable. The color scale maps are identical 
for the C and A run. The velocity vectors for the level are plotted on top of the SSH field. 
It is evident from these maps that the zonal velocity structure has changed between the C 
run and the A run, with three zonal frontal like structures appearing in the A run both at the 
surface and at depth. The temperature structure has also been changed at 160 m showing a 
tongue of higher temperature water at 37°N, corresponding to the higher height field at the 
surface at the same location. The 37.5m layer shows changes also in the temperature 
structure, but the tongue is not as distinct in this layer, although there is a temperature 
increase in the A run as compared to the C run. At both levels, the salinity does not 
appear to have changed as much as the temperature. This is due, in part, to the differences 
in the vertical correlation functions (see Figures 29a,b). It is evident, then, from these 
maps, that not only has SSH changed in this region, but also temperature and salinity, and 
therefore, indirectly, the velocity field. 
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Figure 52. Plan view of model prognostic quantities at level 2 (37.5 meters) in the 
Kuroshio region for the month of December of 1993. The top row (A) are from the 
assimilation run and tlie bottom row are from the C run. The color scale range for each pair 
of quantities is indicated below the panel in the top row. 
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Figure 53. Same as Figure 52, but for level 6 (160 meters). 
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In August of 1992, WOCE investigators (Taft, 1995) collected full depth 
hydrographic data along the 165°E meridian between about 35° and 50°N (PI 3). The A 
run was not begun until after this date, so a comparison is made with the model data from 
August of 1993, with the possibility that the changes between the consecutive Augusts may 
not be that different. Figure 54 (taken from Taft, 1995) shows the vertical structure of 
potential temperature for the WOCE line. Included on the figure are annotations of the 
major frontal structures which were crossed during the sampling of this meridian in the 
Pacific. At around 47°N is the Subarctic Current, between 40° and 42°N is the Subarctic 
Front, and at 35 °N is the location of the Kuroshio Extension. 




Figure 54. From Taft (1995), showing the meridional P13 temperatures (165°E in the 
North Pacfic). 
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Figure 55 shows the vertical temperature structure along this same meridian for the C run 
(Figure 55a) and for the A run (Figure 55b) and for the August Levitus 94 climatology 
(Figure 55c). The heavy lines on the plots indicate the 14*C, 8'C, and 6'C isotherms. 
Note that there are significant changes to the structure in the vertical. It is immediately 
apparent that the 14*C isotherm has moved northward, intersecting the surface at around 
40°N in the A run while in the C run, the isotherm is 5'S at about 35°N. The Levitus field, 
which is being used as the subsurface observational field, also has this isotherm further to 
the north than it is in the C run. The subtropical front at around 40*N shows the isotherms 
closer together in the A run than in the C run. A strip of colder water has been introduced 
at around 43 °N, which looks to be the cause of these compressed isotherms to the south. 
In the north, the 3.5 isotherm is higher in the water column in the A run at 54°N, similar to 
the observations. In addition, the 3-4°C isotherms in the north are more similar to the 
observations in the A run than in the C run. The assimilation process is attempting to 
adjust the subsurface values to those shown in the Levitus field. However, neither model 
run shows the lifting of the deep isotherms south of 37°N. (The observational salinity 
section from PI 3 is not available, therefore it is not shown.) 

b. Gulf Stream Region 

Starting again by showing two levels of plan views of five of the prognostic 
variables of the two model runs (Figures 56: level 2, 37.5 m and Figure 57: level 9, 435 
m), it is easily seen that changes are occurring near the surface as well as at depth. The 
close examination of the path of the Gulf Stream shows that the assimilation run is 
attempting to correct, in temperature as well as in the velocity field, the path of the Gulf 
Stream (Figure 44). A meridional cross section at 70°W is shown in Figure 58a-c for the 
temperature fields from June of 1994 (a) C run, b) A run, and c) Levitus 94). The figure 
shows the northward shift of the 20°C isotherm in the A run. 
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Figure 55. Meridional section along 165°E from a) C run, b) A run, and d) Levitus 94. 
Bold lines indicate the 14°C, 8°C, and 6°C isotherms. 
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Model sections from April of 1994 at 64°W are compared to the April 1985 
Endeavor 129 hydrographic section in Figure 59. Although the years are different (model 
and observations), the comparison is still valid because it is comparing a single month, 
rather than a climatology (e.g. Levitus '94). For this plot, April of 1994 model 
temperatures are shown, with the A run as the bold line and the C run shown by the thin 
line. The in situ data were retrieved from the WOCE hydrographic archive site. The 
observational temperatures, interpolated to 222m (a) and 435m (b), are plotted as asterisks 
to show the difference in the temperatures at these two depths between the model runs and 
the observations. The two models are fairly similar in their representation of the 
temperature fields, especially at the 222m, and both are close to the observed values. 
Where there are differences, near 38°N and at 33°N, the A run is closer to the observed 
temperatures than is simulated in the C run. At a deeper depth of 432m, neither model run 
is simulating correctly the observed temperatures north of 35°N, however the A run 
appears to be adjusting the temperature values towards the observations. This is due to the 
use of a mean field from an earlier model version, and not the true mean of the ocean. For 
the salinities along this section (Figure 60), the assimilation is only slightly modifying the 
values of towards the observed field. 

c. Tropical Region 

A temperature section across 156°E has been extracted from the two runs of 
the model (Figure 61) for two months, July of 1993 (a,b) and six months later in January 
of 1994 (d,e). The corresponding Levitus 94 section is shown in the c and f panels of the 
figure. The sections for the C run are shown in Figure 61a and d, and the sections for the 
A run are in b and e. An examination of the meridional sections shows little change has 
occurred in the region of the western tropical Pacific. 
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Figure 56. Same as Figure 52, but for the Gulf Stream region. 
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Figure 57. Same as Figure 56, but for layer 9 (435 meters). 
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Figure 58. Meridional temperature section along 290°E in the North Atlantic for June of 
1994. Bold lines indicate the 14°C, 8°C, and 6°C isotherms, a) C run temperatures, b) A 
run temperatures, and c) Levitus 94 climatology for June. 
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222m 64W 




Figure 59. Temperatures along section shown in Figure 58 for April of 1994 at a) 222 
meters and b) 432 meters. The bold line are the values from the A run, the thin line are 
from the C run and the astericks represent the observations from the hydrographic cruise 
Endeavor 129 (sampled in April of 1985). 
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Figure 61. Meridional temperature section (156°E) in the tropical Pacific between +/- 8° for 
the months of July 1993 and January 1994. a) and d) are for the C run; b) and e) are the A 
run sections, and c) and f) are the Levitus 94 climatology. 

A closer examination of temperatures across the level of 160 m are shown in Figure 
62. The thin lines represent temperature values for July of 1993, while the bold lines are 
for January 1994. The A run temperatures (solid lines) are cooler as compared to either the 
Levitus data (dotted line) or the C run temperatures (dashed lines) for both months. 



121 



Observations taken from the TAO set of moorings at 156“E between 5°S and 5°N can be 
compared also to the model estimates to determine which run is simulating the observations 
better. Using the gridded maps of the TAO mooring data taken from the PMEL TOGA- 
TAO Web page (www.pmel.noaa.gov/toga-tao), approximate monthly mean temperature 
observations plus the associated anomaly are also plotted on this figure. The July 1993 
mean temperature at the given latitudes and at 160 m are plotted as circles and January 1994 
temperatures are plotted as asterisks. The model July 1993 temperatures are still too warm 
as compared to the observations (compare the thin lines to the circle values), but the A run 
has reduced the difference by approximately a degree. The A run at 5“N is accurately 
representing the ocean temperature, once again reducing the temperature from what is seen 
in the section from the C run. At 0°N, the A run has reduced the temperature, but not to the 
mean value seen in the observations. And in the south, the two runs differ by only a few 
tenths of a degree, with the C run closer to the observations than the A run. It is noted that 
the monthly mean model temperatures from both model runs are within the range of the 
temperatures seen in the observations over these two months. The Levitus climatology 
shows warmer values than the observed or modeled values. For the period 93/94, the 
observed values are cooler than the average values' as seen in time/latitude plots of TAO 
data covering six years. Thus the model is making the adjustment to the temperature fields 
in the correct direction. 

d) Southern Ocean 

WOCE line All {Saunders and King, 1995), along 45°S between 60°W 
and 20°E in the Atlantic section, was sampled December, 1992 - February of 1993. It is in 
the earlier part of the assimilation experiment, but it can still provide some guide as to how 
realistic either of the model runs are. Zonal sections across 45°S are shown in Figure 
63a,b for January of 1993 and in Figure 64a,b for January of 1994 for run C and run A, 
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respectively. The January 1993 plots, in the early part of the experiment, for the C (Figure 
63a) run and for the A run (Figure 63b) are similar, with the major difference at 50°W, 
where the large plume like structure of warm water has increased in the A run. These can 
be compared to Figure 63c, the A1 1 observational temperatures sampled around January of 




Figure 62. Tempertures at 160 meters for the section shown in Figure 61. The solid lines 
(bold = January, thin = July) are the A run temperatures, the dashed lines are from the C 
run and the dotted lines are the Levitus values. The astericks denote the January 1994 TAO 
observations and the circles, July 1993 TAO observations at the same depth and location 
with their associated variance. 



1993. These plots are directly comparable. In the observations, the plume of warm water 
is located at 56°W rather than at 50°W as seen in the model. Another obvious difference is 
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Figure 63. Zonal sections of temperatures along 45°S in the Atlantic Ocean (between 60'W 
and 0°E) for January of 1993 for a) the C run, b) the A run, and c) the A1 1 observed 
temperatures. 
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Figure 64. Zonal sections of temperatures along 45°S in the Atlantic Ocean (between 60°W 
and 0°E) for January of 1994 for a) the C run, b) the A run, and c) the variance map of the 
yearly Levitus climatology. 
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that the water at depth in the observations (Figure 56c) is cooler than either of the model 
runs. The observations also show fewer of the wave-like temperature structures visible in 
the model simulations. By January of 1994, the A run (Figure 64b) has developed many 
deep penetrating features, while the C run has remained approximately the same. Thus, 
this suggests that the wind field is quite similar between the two years because the C run is 
similar in both years, and the changes seen in the A run are purely due to the assimilation 
process. 

In this region, the assimilation run is very unrealistic and the source of the model 
diverging away observations is the variance field seen in the monthly Levitus fields along 
this line (Figure 64c). It can be seen that between 50°W and 55°W, the variability of the 
Levitus data is relatively high. This, along with the corresponding functional, which 
shows a highly barotropic forcing function for this region down to a depth of 1000 meters 
or so, is responsible for the modifications of the model fields. Not much observational 
data has been collected in this region and the Levitus climatology may not be realistically 
representing the monthly tenperature fields. 

4. Heat Content and Transport 

The heat content of the ocean is used as an indirect measure of how the ocean is 
transporting water, and thus heat around. Accurate or close representation of the heat 
content and transport in the ocean is necessary prior to coupling ocean models with 
atmospheric models to study climate change. Figure 65 shows the northward meridional 
heat transport in the model computed from the 2 year mean (1993/1994) for the C run 
(dashed line) and the A run (solid line). The global northward heat transport is shown in 
Figure 65a, with the basin integrated transports shown in b), for the Atlantic and c) for the 
Indo-Pacific region. The two runs have produced similar transport values. In the Atlantic, 
the two runs are virtually identical. Most of the difference in the values shown in the global 
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view is due to the differences seen in the Indo-Pacific basin. Here, run A has reduced the 
transports in both the northern and southern hemisphere. 

There have been many studies which have attempted to estimate heat transport from 
the observational data. The heat transport values from a recent and comprehensive study 
by Macdonald and Wunsch (1996) are plotted on Figure 65, for the globe and both the 
basins as circles with an error bar. The model (both the A and C runs) is reproducing 
accurately the transports that are observed in the Atlantic. In the Indo-Pacific, Figure 65c, 
the A run heat transports are closer to the observed values in the northern hemisphere. In 
the southern hemisphere, the C run is closer, although the modeled values still deviate from 
the observations. Macdonald and Wunsch have also provided some estimates of heat flux 
across several lines in the Southern Pacific and Indian Oceans. These data (Table 5) show 
that it is the Indian Ocean (rather than the Pacific) in the model which is producing most of 
the difference between the observed transports and the modeled transports and between the 
2 modeled transports at 32°S. 



Section 


Macdonald 
and Wunsch 


C run 


A run 


Pacific 28°S 


-0.04 ±0.30 


0.26 


0.38 


Pacific 43°S 


0.26±0.28 


0.33 


0.37 


Indian 18°S 


-1.45±0.30 


-0.37 


-0.31 


Indian 32°S 


-1.3010.28 


0.87 


1.13 



Table 5. Heat flux in petawatts across specific zonal sections in the Indo-Pacific Oceans. 

The amount of heat at a location is examined at various locations around the globe. 
Plots of the monthly cycle of the integrated heat content of the upper layers of the model (to 
310m) with respect to the surface temperature are shown in Figures 66 (A run), 67 (C run), 
and 68 (Levitus 94). Monthly value of heat content are plotted separately for 1993 and 
1994 Figures 66 and 67, with the bold line representing 1993. A realistic model should 



127 



show a hysteresis effect in the heat content. This means that the heat content of the ocean 
for a particular month should retain some of the heat of the prior month and should not be 
directly related to the changes seen in the surface temperatures. In other words, as the 
ocean cools during the winter months, some of the heat gained during the warmer months 
should still be retained in the water column. This is shown in the northern hemisphere by 
the months of September through February having higher heat content than in the months 
of August to March, even though the surface temperatures are similar. (The loop starts 
with December of the year before, denoted by the heavy dot.) The assimilation process 
(Figure 66) corrects the model's heat content to some extent, especially in 1994. In effect, 
the assimilation is forcing the model's temperatures to be more representative of the mixed 
layer seen in observations (somewhat represented by Figure 68, the Levitus data). If a 
mixed layer formulation is included in the model, this deficit of the model can be overcome 
without the need for assimilation (pers. comm. Semtner). 

A closer look at one location, 35°N, 30°W (Figure 69) shows the differences in the 
monthly structure of the water column between the two runs of the model and between the 
two years of the model runs. In both years, beginning in February, the column is cooler in 
the A run than in the C run, with the 16‘’C isotherm surfacing in March of 1994. There is 
little difference in the surface waters (top 70m) during the second half of the year. In the 
winter months, the 16°C isotherm uplifted, showing a deeper "mixed" layer, especially in 
1994. 
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Figure 65. Meridional northward heat transport for the 2 year mean (93/94) for a) the 
globe, b) Atlantic, and c) Indo-Pacific. The bold line is the heat transport computed from 
the A run and the dashed line is the transport from the C run. The circles denote the 
Macdonald and Wunsch (1996) values computed from observations. 



129 



-34.37N219.4ES.Pac. 



-34.37N 319.4E S. Atl. 





35.03N 330.2E N. Ati. 




surface temperature degree C 





53.46N319.4EN. Atl. 




surface temperature degree C 



Figure 66. Hysteresis curves for heat content in run A at six locations; a) 34.37°S 219.4°E 
in the South Pacific, b) 34.37°S 319.4°E in the South Atlantic, c) 16.38°S 59.4° E in the 
Indian, d) 35.03°N 179.4°E in the North Pacific, e) 35.03°N 330.2°E in the North 
Atlantic, and f) 53.46°N 319.4°E in the North Atlantic. The bold curves are for 1994 and 
the thin lines are for 1993. The circles denote the starting location and x marks the months. 
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Figure 67. Same as Figure 66, but for the C run. 
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Figure 68. Same as Figure 66, but for the Levitus climatology. 
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Figure 69. Time/depth plot of temperature at 35°N, 30°W for the a) 1993 A run, b) the 
1994 A run, c) the 1993 C run, and d) the 1994 C run. 

5. Thermohaline Circulation 



The most easily detected changes resulting from the application of the assimilation 
process occur at or near the surface because the observational field that is related to the 
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surface quantity, SSH. Indirectly, what occurs at the surface will, in time, effect the deep 
circulation because of the changes made to the density of the water. This thermohaline 
circulation, relating to modifications to the temperature and salinity field, are seen in nature 
when, for example, there is change in the ice melt at high latitudes or changes in the 
evaporation-precipitation rates over an ocean. The assimilation experiment not only affects 
surface properties, but changes temperature and salinity at depth via the correlations of 
SSH with the vertical structure of temperature and salinity, in a more accelerated maimer 
than occurs in nature. Thus, these changes at depth, should also change the thermohaline 
circulation of the model. The thermohaline structure will only be discussed in terms of the 
changes seen in the two model runs, since complete, global and detailed observational data 
is not available. Global WOCE data, covering about the same period of the model runs, 
should soon become available to make a complete comparison of the observational 1990s 
thermohaline circulation. 

The meridional overturning streamfunction is one method to examine the 
differences in the model runs of the thermohaline circulation. Plots of meridional 
overturning show how water circulates between the north and south basins of an ocean. 
The streamlines can be used to trace the flow from where it enters at the surface down to 
the subsurface layers of the ocean, as is described by the ocean "conveyor belt " system. 
The equation for meridional overturning is: 



W(?t,0 = 



J v(<p,X,z') a cos A dz'd<p 

z(<t>,X,0 



4.6, 



where X and <}) are latitude and longitude, respectively, and a is equal to the radius of the 
earth. ^ is the depth of the surface being considered, either a density surface or the 
traditional level surface. The model meridional velocity is first integrated vertically to the 
specified surface (either level or density), then zonally integrated for the globe or a basin. 

The Indo-Pacific, Atlantic, and global meridional overturning streamfunction 
defined on level surfaces is shown in Figure 70a-c, with the A run shown in the top panel 
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of each figure and the C run on the bottom for the mean of 1994. From the global 
meridional integrated sections, one can see the Antarctic cell (0-2000m, 65°S), the Deacon 
cell (between 60°S and 40°S and 500-3000m), the deep cell of northward following water, 
the shallow subtropical gyre circulations in both the north and the south and the mid-depth 
flow which represents the major transfer of the thermohaline transport of water from one 
basin to another. 

There are several immediate differences seen between the model runs. In the Indo- 
Pacific region (Figure 70a), the major differences are in the flow of the intermediate waters 
at 20°S, defined by the shaded areas of the southern hemisphere. The assimilation run (A) 
shows strong northward flow (average of about 6 Sv), ventilating the interior of the Indian 
and/or the Pacific ocean to about 10°S, and then returning in the deeper layers. Most of the 
northward flux of water is in the Indian ocean as the suggested from the heat flux 
calculations (Table 5) shown in the previous section. This northward flow is nonexistent 
in the control run of the model. In the North Pacific, the mid-depth southward flow 
centered at 25 °N in the C run is much reduced, or eliminated in the A run. The -2 Sv 
contour above 600 meters is much stronger in the A run than in the C run. In fact, it is 
similar to the streamlines calculated for the 0.5° resolution global ocean model run of the 
Semtem/Chervin model {Semtner and Chervin, 1992, Figure 9). The 0.5° model was 
constrained to Levitus data whereas the 1/4° resolution run was not (i.e., the run C). The 
assimilation, in some fashion, also constrains the model to a Levitus data set, so it is not 
surprising that the streamlines of the A run resemble this earlier 0.5° model run. 

The Atlantic region (Figure 70b) shows similar circulation patterns above 600m. 
In the deep ocean, too much forcing was applied at 35°N, where the flow looks somewhat 
unrealistic in the A run. The deep northward flow visible in the streamlines from the A run 
(denoted by the shaded region) is the Antarctic Bottom water. This water is not seen in the 
C run of the model (or is much smaller, below 2Sv.) Standard oceanography texts, such 
as Pickard and Emery (1986) show a similar deep flow. Figure 70c. showing the global 
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view of the overturning streamlines is included for completeness and to define where the 
streamlines for the deep flows in the southern hemisphere begin. 

Following Hirst, et al. (1996) and Dods and Webb (1994), overturning flow can 
also be shown along surfaces of constant density. For this comparison, the 03 surface has 
been chosen because it is very similar to the representation of streamlines on neutral density 
surfaces (Hirst, et al.). Figures 71a-c show Indo-Pacific, Atlantic and global meridional 
overturning on O3 surfaces calculated from both the A (top panels) and C (bottom panels) 
runs of the model. Viewing streamlines on density surfaces shows the flow in a more 
realistic manner. The major recirculation cells again can be seen. For example, the Atlantic 
circulation associated with the transfer of water from the northern to southern hemisphere 
("the conveyor belt") can be seen in Figure 71b, on the density surface of about 41. to 
41.5. The panels for the Indo-Pacific region show similar features to the plots on the level 
surfaces. The A run has somewhat weaker circulation around the northern basin (shaded 
area) than the C run, and the deep flow is stronger in the A run than in the C run. The A 
run has removed what appears to be a narrow and deep penetrating feature (6Sv contour) 
located at 12°N. The deep northward flow of the Antarctic bottom water is visible, 
flowing northward to about 20°N, not quite as far north, as seen when plotted on level 
surfaces (Figure 70b). 

Figure 71c, the global overturning streamfunction shows the intense equatorial 
features which are weaker in the A run, than in the C run. The Deacon cell, is denoted by a 
very small circulation cell, centered at about 50°S and on a density surface of 41 in Figure 
71c. 



6. Summary 



This section of Chapter IV has shown, with both qualitative and quantitative 
comparisons made to observational data from altimeters, hydrographic sections, and tide- 
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gauges, the assimilation process is adjusting the model's prognostic variables towards 
observed values. These adjustments can be seen both in the surface fields and subsurface 
in the model's variables which have been directly nudged (SSH, temperature, salinity) and 




Latitude 



Figure 70a) Meridional overturning streamfimction defined on level surfaces for the Indo- 
Pacific region with the top panel showing the A run and the bottom panel, the C run. 
Contours are every 4 Sv. Shaded areas define regions referenced to in the text. 
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Figure 70b) Same as Figure 70a, but for the Atlantic. 
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Figure 70c) Same as Figure 70c, but for the globe. 
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Figure 71a) Meridional overturning streamfunction defined on 03 surfaces for the Indo- 
Pacific region with the top panel showing the A run and the bottom panel, the C run. 
Contours are in SV 10e6 kgm"^ and are shown by increment of 4 SV. 
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Figure 71b) Same as Figure 71a, but for the Atlantic. 
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Figure 71c) Same as Figure 71a, but for the globe. 

in the variables which are indirectly modified (velocity). Significant changes can also be 
seen in the diagnostic quantities such as meridional overturning and heat content. The 
assimilation process fails in those regions (mainly the Southern Ocean) where there is a 
lack of in situ data to correctly adjust the model's subsurface fields. These areas which 
under perform can be better simulated once the observational fields are improved. 
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E. CONCLUSIONS 



The assimilation process for the assimilation SSH and associated (via correlation 
functionals) temperatures and salinities has been described. The model's prognostic 
variables from the same period of time in the model run are statistically different between 
the assimilation run (A) and the control run (C) as shown with the use of the Student-t test 
and the F-test, and with the comparison of the mean and variance of SSH, temperature, and 
salinity. These differences are both at the surface and subsurface, showing that the 
assimilation process is making changes in both, though to a different degree. It has also 
been shown that the method is dynamically compatible with the numerics of the primitive 
equation ocean model. And finally, the synthesized ocean variables, in most of the regions 
with high variance, have been adjusted towards values seen in observations. 
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V. IMPACT OF DIFFERENT SAMPLING PATTERNS IN 

OBSERVATIONS 

Altimeter data from two different satellites have been used as observational SSH 
information for the second assimilation experiment. For 1993, both ERS-1 and T/P data 
were used, while in 1994, only T/P data was available. Each satellite has a different 
sampling period, ERS-1 at a 35 day repeat and T/P, with a 10 day repeat cycle. This 
provides an opportunity to quantify the advantage of having multiple satellites flying with 
different orbits over having just one sateUite. 

Statistics are used to quantify the difference in the two years of the assimilation run 
to establish which year was more successful at incorporating the observational SSH signal 
into the model. Statistics have been computed to determine 1) which year is closer in its 
representation of the assimilated SSH field and 2) which year has increased the variance of 
the SSH the most. The second question is asked because the model is known to be under 
representing the SSH variability when it is mn without any assimilation. 

For the first question (which is closer to the altimeter field it assimilated), the mean 
square difference (MSD) over all points is calculated between the monthly mean SSH of 
run A and the monthly mean of the altimeter data, as defined in equation 5.1. 



MSD-1 = [<M> - <D>]2 j j 5.1, 

where <M> is the monthly mean SSH of the model and <D> is the monthly mean SSH of 
the assimilated altimeter data at a location i,j. To compare the two years, 1993 and 1994, 
the MSD values calculated for each year are normalized by the variance of <D> over all i,j, 
to account for differences in the number of points at which observations were assimilated at 
for the two years. Equation 5.1, as used for this comparison is defined as: 
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MSD-1 = [<M> - <D>]2 i j / Var(<D>) 



5.2. 



The two numbers from each year are compared, with the smaller number indicating which 
year is closer to the observations that were assimilated. Table 6 shows computed MSD 
values for the two years. For all the months, the values are smaller in 1993 than in 1994, 
suggesting that the increase in the number of observed locations (ERS-1 samples at a finer 
spatial resolution) is beneficial in the assimilation process. 



Month 


MSD-1 93 


MSD-1 94 


JAN 


1.77 


2.16 


FEB 


1.82 


2.16 


MAR 


1.77 


2.15 


APR 


1.82 


2.16 


MAY 


1.80 


2.17 


JUN 


1.83 


2.20 


JUL 


1.82 


2.22 


AUG 


1.84 


2.23 


SEP 


1.86 


2.25 


OCT 


1.84 


2.24 


NOV 


1.86 


2.23 


DEC 


1.99 


2.18 



Table 6: Monthly normalized mean square difference between run A SSH and Observed 
SSH. 



Similarly, to examine whether or not the model's variance has increased, instead of 
computing the mean square difference between the assimilated run and the observations, 
the MSD is computed between monthly variances of run A and run C, so that equation 5.2 
is rewritten as: 



MSD-2 = [<M> - <C>]2 i j / Var(<C>) 



5.3, 



where <M> is now the monthly variance at a grid point for mn A and <C> is the same for 
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run C. The higher the value of MSD-2, the higher the difference in the variances between 
the two model runs, A and C. Table 7 shows the computed values. Eight of the twelve 
months have values higher in 1993 than in 1994, with two months (July and April) fairly 
similar and two months where the C run has the higher increase in its variance (January and 
May). It is not surprising that January of 1994 is higher than January of 1993, since 
January of 1993 is fairly early in the model run. 



Month 


MSD 93 


MSD 94 


JAN 


0.43 


0.49 


FEB 


0.94 


0.58 


MAR 


0.63 


0.45 


APR 


1.09 


1.16 


MAY 


1.14 


1.47 


JUN 


0.96 


0.76 


JUL 


1.42 


1.48 


AUG 


1.07 


0.82 


SEP 


1.25 


0.66 


OCT 


0.70 


0.50 


NOV 


0.66 


0.45 


DEC 


0.87 


0.53 



Table 7: Monthly normalized mean square difference between mn A SSH and run C SSH. 

To summarize, the above quantitative evaluation suggests that it is useful to use 
both ERS-1 data with its higher spatial sampling and T/P with its higher temporal sampling 
as the observational data set when ever possible. This has been determined by the lower 
mean square difference between the model SSH fields and the observations for 1993 (verse 
1994). In addition, mn A of the model also has a higher increase in variance in 1993, than 
in 1994, suggesting that it is the finer spatially sampled ERS-1 data, which includes 
mesoscale features, influencing this variance in the model. 
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VI. POSSIBLE MODEL IMPROVEMENTS 



As has been discussed in the previous chapters, when the ocean model is run 
without the assimilation of any observations, the model simulates the ocean adequately. 
The model is able to do this because its resolution is high enough to resolve narrow and 
fast currents, such as the western boundary currents (Gulf Stream and Kuroshio). The 
model has a free surface and well defined bathymetry which has not been smoothed and 
which influences the path of the currents. The model also uses monthly varying heat fluxes 
producing a seasonal signal in the model's SSH (Stammer, et al, 1996). Most 
assimilation experiments previously have not used a control run which include these 
aspects, thus producing a control run which is much less accurate. For example, Oschlies 
and Willibrand [1996] use monthly mean wind stresses computed from climatological data 
sets for their assimilation experiment in the Atlantic ocean and Carton, et al. [1996] use a 
relatively low resolution PE model (0.5° x 2.5° grid in the tropics) and weekly 1000 mbar 
National Centers for Environmental Prediction (NCEP) wind fields rather than the higher 
temporal varying 3-day 10 m ECMWF winds, which are used in these experiments. 
Therefore, because this experiment is using more accurate forcing (compared to other 
assimilation experiments) and improved numerics (i.e. the inclusion of the free surface as a 
prognostic variable), the modifications made by the assimilation to adjust the model 
towards reality are relatively small. This chapter looks at where changes can be made in 
the model's numerical formulation, forcing fields, or initial conditions to better simulate the 
ocean circulation without the need for assimilation. 

A. IMPROVEMENT IN THE MODEL'S MIXED LAYER 
FORMULATION 

The assimilation process adjusts the model's temperature field and in the process. 
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accounts for the temperature changes which would be seen if a numerical code to handle 
mixed layer dynamics were included in the model's formulation. This is the clearest of the 
assimilation effects on the model to understand. The bold line on Figure 72 shows the 
histerisis curve at various points from 1994 of run A, reproduced from Figure 66. The thin 
lines, on the same figure, represent the curve as calculated from an initial experiment which 
includes mixed layer dynamics in the model formulation (pers. comm., Semtner, 1997). 
The mixed layer formulation is the Krause-Tumer {Krause and Turner, 1967) type. The 
two model runs do not cover the same period. The A run is for 1994 and the mixed layer 
run covers a time period from Sept, of 1992 through August of 1993. There are 
similarities between the curves, especially at 35“S, 219°E in the South Pacific and 53.5°N, 
3 19°E in the North Atlantic, where both differ from run C (see Figure 66). 

The differences in the temporal change of the temperature fields over twelve months 
are shown in Figure 73 with panels showing the temperature with depth (to 300 m) for 
each of the six locations shown in Figure 72. As an example, in the North Atlantic at 
53.5°N, 319°E (last row of panels on Figure 73), both the A run and the mixed layer run 
show increased temperatures at depth starting in about March over the temperatures seen in 
the C run. In the South Pacific (first row of panels on Figure 73), the deeper layers (below 
100 m) cool between January and December for the A run and the mixed layer test, while 
the C run temperatures are fairly constant. 

These figures, thus, suggest that assimilation can compensate for a model lacking in 
mixed layer dynamics by adjusting the temperatures in the upper layers to account for 
seasonal changes. There are inaccuracies in both model runs. Mixed layer run was an 
initial test of the numeric code and still needs improvements. At this point, it can not be 
said that mixed layer code will completely correct inaccuracies in the model's top layers and 
eliminate the use of assimilation to correct this deficit. 
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Figure 72. Heat Content for 1994 of the A run (bold line) and Mixed layer run, 1992/3 (thin 
line) at six locations in various ocean basins. 
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Figure 73. Temporal temperature change over 12 months at 6 locations for the A run, the C 
run, and mixed layer test. Contours are 1 °C. 
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B. IMPROVEMENT OF THE MODEL'S INITIAL CONDITIONS 



Adjustments to the model's initial conditions can be made by examining the changes 
made to the mean over the two year period, 1 993-1994. Figure 32a (Chapter IV) showed 
the difference in the mean field between the two model mns (A and C). The region to the 
east of Japan has a large change in its mean field at the surface (Figure 32a) and also in the 
temperature field at 160m (Figure 33a). If a meridional section of temperature and salinity 
is plotted along 165°E between 30° and 45°N, from the 2 year mean fields of the A run and 
C run, along with the Levitus 94 climatological data (Figure 74), run A qualitatively 
appears to be adjusting towards the Levitus data which has been used in the creation of the 
vertical functionals which relate SSH to temperature. The 14°C and 15°C isotherms in the 
A run are sloping upwards (from left to right ) similar that shown in the Levitus data, while 
the C run has a very sharp vertical front near 34°N. The RMS difference between the 
temperature fields are: A mn and Levitus, 1.80°C (to 200 m); C mn and Levitus, 1.96°C. 
For the salinity fields, the numbers are; between the A mn and Levitus, 0.1 l%o; between C 
mn and Levitus, 0.07%o. The average T/S diagram of the upper layers also shows and 
adjustment towards the Levitus data (Figure 75). 

Similarly, a section can be looked at in the Indian Ocean, along 40°S between 30°E 
and 120°E (Figure 76). The western side shows large differences between mn A (Figure 
76b) and mn C (Figure 76c), in the region of the Aghulas Current. 

It was noted in Table 5, that across 32°S in the Indian Ocean, the heat flux is 
underestimated by both model mns, with a larger error seen in mn A. An examination of 
the sections across 32°S (not shown) and 40°S, indicates that the difference between the 
two model runs is on the west side. The assimilation pushes the model towards the Levitus 
94 mean (used as the observational field for the vertical functionals), which results in the 
weakening of the current southward and the heat flux balance increases northward. 
(Levitus shows little or no current, due to the low resolution of the climatology.) This is 
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Figure 74. Meridional sections of temperature (left panels) and salinity (right panels) along 
165°E in the North Pacific for a) Levitus 94 data, b) A run, and c) C run. Contours are 1°C 
for temperature and 0.1 S for salinity. 
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shown in Figure 77, which shows the T/S diagrams for the full zonal section (a), in the 
Aghulas current area (b) between 33.5°E and 37.5°E, and mid way along the section (c) 
between 65.5°E and 69.5”E. The RMS difference between the Levitus mean and the run A 
is 1.37°C for temperature and O.l6%o for salinity and for run C, 1.4rC and 0.17%o, 
meaning that A is closer to Levitus than run C. The mean flow of either the model or that 
derived from the Levitus 94 data, is not correct, and thus limits the assimilation process in 
correcting for these errors in the model simulations. If given a different mean or initial 
condition, with associated confidence levels, which is more refined and accurate, the 
assimilation process can be used to adjust the model's prognostic variables towards this 
initial or mean state. The vertical structure which relates SSH to T and S does not need to 
be defined at high resolution over the globe, but should be at places which have known 
narrow currents. 




Figure 75. T/S diagram for meridional section shown in Figure 74. Bold line is for Levitus 
94 data, thin line represents the A run values, and the dotted line represents the C run. 
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Figure 76. Same as Figure 74, except for 40°S in the Indian Ocean. 



156 




Figure 77. T/S diagram for zonal section shown in Figure 76. Bold line is for Levitus 94 
data, thin line represents the A run values, and the dotted line represents the C mn. 

These examples are encouraging, given that the observational SSH field includes a 
mean a similar to the mean of the control run and uses a coarse monthly climatology to the 
surface SSH to subsurface temperature and salinity fields. More accurate adjustment 
should occur if the vertical relationships are computed using data at higher spatial 
resolution in areas where narrow currents occur. In addition, repeat WOCE hydrographic 
data can completely replace the empirical vertical T/S profiles with that of real data. 

C. IMPROVEMENT IN THE TEMPORAL VARIABILITY AT A 
LOCATION 

The changes in the temporal variability of the ocean at a location has been 
previously examined when the time series of SSH from the model runs were compared to 
similar series obtained from tide-gauges. Another way to look at changes related mostly to 
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mesoscale processes is to look at a zonal section of SSH in time. Figure 78 shows such a 
section at 34°N in the Atlantic between 30°W and 0°, a region with a medium weighted 
nudging. The temporal mean SSH has been removed, along with an aimual harmonic. 
Run A (top plot) shows stronger westward propagating waves in the eastern half of the 
domain than is visible in the plot for run C (lower plot). These westward propagating 
Rossby waves have been seen in altimeter data beginning with Geosat data {Tokmakian and 
Challenor, 1993), in addition to being observed in sea surface temperature maps (Cipollini, 
et al. 1997). Plotted in Figure 79, are maps of SSH for the region which can be used to 
compare to the observations gathered in July of 1993 (day 201) and September 1993 (day 
250) by the French Semaphore experiment {Hernandez and Le Traon, 1995). A spatial 
mean over the area has been removed to more easily compare to the Hernandez and Le 
Traon's Figure 6. The A run (Figures 79a and c) show a depression in SSH centered at 
25°W,35°N in July (day 201) which has moved to the west by day 250. The SSH 
structure is more closely resembling the hydrographic surface height than the SSH maps 
from ERS-1 and T/P shown by Hernandez and Le Traon. This is due to the having a mean 
SSH from a previous model run added to the SSH altimeter anomalies used in the 
assimilation process, rather than a true mean. It may also be due to a weak nudging 
coefficient. 

This example shows that mesoscale features can be enhanced and introduced into 
the model. Other examples can be found in the flow of the Antarctic Circumpolar Current. 
Warm and cold core eddies of the Gulf Stream and the Aghulas are not as easily identifiable 
because the assimilation process needs to correct the mean path of the Gulf Stream (or 
Aghulas), in addition to correctly placing an eddy feature or removing an eddy. 
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Figure 78. Time/longitude plot of SSH along 34°N in the North Atlantic in the region of the 
Azores Current, a) A run and b) C run. 
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Figure 79. Plan maps of SSH for the North Atlantic in the region of the Azores current, a) 
Julian Day 201, 1993 for the A run, b) Julian Day 201, 1993 for the C run, c) Julian Day 
250, 1993 for the A run, and d) Julian Day 250, 1993 for the C run. 



D. SUMMARY 



The assimilation has changed the model, producing an ocean which includes some 
of the properties associated with a model which has included a mixed layer in its numerics. 
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Additionally, the assimilation has changed the model's mean state, which suggests that 
given a more accurate initial state, the model will reproduce observed and derived 
quantities, such as heat transport, more accurately. Each of these model "corrections" 
produced a model that is dynamically stable and compatible with the control run. This 
suggests that once a mixed layer is included in the numerics and an improved initial 
condition at high spatial resolution is incorporated, the assimilation process will not change 
these aspects of the model's output. 

The shorter term variations in the ocean may not be as easily reproduced without 
having included in the forcing the event that caused the variability. For example, will a 
more accurate wind stress field in the region of the Azores Current reproduce the westward 
propagating Rossby waves or does the model need the assimilation to correctly "time" the 
events? It may never be possible to have a surface forcing field that is accurate enough to 
reproduce these events, thus if a model is being used to forecast the location of a front, or 
an eddy, assimilation is necessary. But if the model is used for climate purposes, if the 
events are statistically accounted for over a period of several years, than assimilation is 
probably not warranted. 

In short, this experiment of using observational SSH anomaly measurements as the 
driver for the assimilation fields shows that a model can be modified to account for 
different initial conditions, or not having a mixed layer, and can introduce or enhance 
features that are not seen in the model run without assimilation. 
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VII. CONCLUSIONS 



The research presented in this thesis on the evaluation of two assimilation 
experiments has shown that it is feasible to assimilate a large quantity of observations into 
an eddy-resolving global ocean model. 

First, it has been shown, using a twin experiment, that model diagnostic fields can 
be adjusted directly (SSH, temperature and salinity) and indirectly (velocity) using SSH 
measurements and empirically defined functions which relate SSH to subsurface 
temperature and salinity. After the assimilation of SSH using a simple nudging process, 
these prognostic fields of the twin run resemble the control fields rather than the initial 
fields. 

The analysis of the second experiment (using SSH anomaly fields sampled with the 
ERS-1 and the TOPEX/POSEIDON altimeters) shows that assimilation using the nudging 
method is dynamically consistent with the control run of the model. This experiment has 
shown that the simulation has adjusted the model towards observed values of the global 
ocean. However, the subsurface values are a climatology and do not reflect the true ocean 
variability in time or space. Thus, the assimilation is not successful at recreating fine 
narrow structures, such as the Aghulas Current, accurately. The assimilation process has 
changed diagnostic quantities, and integrated quantities, indicating that assimilation can 
account for some processes, such as the seasonal mixed layer temperature variations, 
which may not be included the numerics of a model. Other analyses, for example the 
meridional overturning, show that assimilation modified the deep abyssal flow, as well as 
the near surface flow. 

At the mesoscale, instabilities have been introduced where observations show that 
they exist, such as the Azores current. Nudging also has attempted to modify other 
mesoscale flow, such as incorrectly placed eddies and features in the region of the Gulf 
Stream. This is in addition to modifying the mean ocean circulation field. 
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This research has shown that a simple, straight forward method of assimilating 
observational data into a high-resolution global ocean model is dynamically consistent 
because the model's own dynamics correct for any dynamical inconsistencies in the 
observations. This is a reasonable method to use and can be further enhanced by replacing 
the Levims 94 monthly climatology with regionally collected observations, especially those 
which have an associated variance. Such data are now, or soon to be, available in the form 
of the WOCE repeat hydrographic sections. Other WOCE data can be used to enhance the 
climatological field that has used data accumulated over many years. Subsurface data can 
also be assimilated in a direct manner, without relating it to the SSH as has been done for 
these experiments. 
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APPENDIX 



The mean SSH field (<SSH_observed>) which is added to the altimeter SSH 
anomalies effects part of the change in the resulting yearly mean SSH field. In the 
assimilation run (run A), a mean different from the control run's mean SSH field is used to 
show what may occur if a mean computed from observations (assumed to be different from 
the control run's mean) where to be used. To understand the effect of this mean SSH field 
on the simulation, a short (approximately 1 year) test (referred to as the E run) was 
conducted using the 2 year SSH mean (93/94) from the C run instead of the 3-year mean 
(87-89) from the POCM_4A run (see Figure 28). Figure 80 shows the difference between 
the two mean fields of Figure 28. Some of the differences are due to the different forcing 
used, with the POCM_4A run using the ECMWF lOOOmb wind stress fields, known to be 
overestimated, while the POCM_4B run (the control run) used the more accurate 10 m 
wind fields. 



POCM_4B 93/94 - POCM_4A 87-89 > 10 cm 




Longitude 



Figure 80. Difference in the mean SSH field for 3 years (87-89) of the POCM_4A run and 
the POCM_4B run (93/94). Red indicates POCM_4A is larger (-10cm) and blue indicates 
where POCM_4B is larger (10cm). (These are not absolute differences, e.g.. if 4A = -200 
cm and 4B = -190 cm, then 4B - 4A = 10cm. The color scheme is being used to show just 
that there are differences.) 
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The next two figures show the effect each of the <SSH_observed> fields has on the 
simulation. Figure 81a and b show the differences in the simulated 1993 SSH mean field 
as compared to the control run (run C) for the E run (<SSH_observed> = POCM_4B, 93- 
94) and for the A run (<SSH_observed> = POCM_4A, 87-89). When the 
<SSH_observed> is the mean of the control run, the differences between the control run 
and the assimilation are reduced, especially in the area of the Kuroshio Extension (Figure 
81a). The mean SSH of the assimilation runs were still modified, though, independent of 
which mean SSH field was added to the altimeter SSH anomalies as shown in Figure 81c. 
Here, the difference between the assimilation runs shows only small differences outside of 
the Kuroshio area (such as in the Gulf Stream region) and no difference between the two 
assimilation runs in the Indian Ocean. Both assimilation runs show changes in the mean 
SSH field occurring in the Indian Ocean even though there is little if no difference in the 
<SSH_observed> field (see Figure 80). 

In the plot of the mean 1993 temperature field at 160m (level 6 of the model), the 
<SSH_observed> has mostly contributed to the change seen in the Kuroshio region of the 
A run (Figure 82c). As with the SSH fields, the mean temperature field at depth is 
modified independent of which <SSH_observed> was used, such as in the Southern Ocean 
between 20° and 70°E, 45°S (Figures 82a,b). 

In summary, both the <SSH_observed> and the SSH anomalies combined with 
the associated vertical adjustments, modify the mean fields of temperature and SSH when 
assimilating observations in the global model. In some cases, such as in the Kuroshio, the 
difference in the mean fields can be directly related to the <SSH_observed> used. Also, 
the simulated mean SSH does not always change, when there is a change in the mean 
temperature field does, as is seen in the temperature field of the tropical ocean in Figure 33a 
verse what is shown in the mean SSH field (Figure 32a). 
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Figure 81. Differences between pairs of simulated mean SSH fields, a) between the C run 
and the E run, b) between the C run and the A run, c) between the two assimilation runs, A 
and E. Red indicates -10cm or smaller {i.e. , more negative) and blue indicates 10 cm or 
larger. 
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Figure 82. Differences between pairs of simulated mean temperature fields at 160m, a) 
between the C run and the E run, b) between the C run and the A run, c) between the two 
assimilation runs, A and E. Red indicates -TC and smaller (i.e . , more negative) and blue 
indicates 1°C and larger. 
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